This R Markdown document describes the analyses performed for the manuscript entitled “Environmental pollution correlates with parasite infection across a riverine landscape” by Io S. Deflem, Seppe Marchand, Federico C.F. Calboli, Joost A.M. Raeymaekers, Filip A.M. Volckaert and Pascal I. Hablützel.
The analyses were run in R 4.2.2
Up to thirty 0+ three-spined sticklebacks were sampled at 37 locations in the Dijle and Demer basins in Belgium during autumn 2016 under a permit of the Flemish Agency Nature and Forest. Both basins together cover a continuous surface area of 3,624 km² with the furthest two sampling sites being located 117 km apart (distance measured along rivers). All locations included small and relatively slow flowing streams (drop off from highest to lowest point is 18 m) covering a wide range of ecological, hydromorphological, and physico-chemical characteristics. Fish were caught using a dip net.
# Empty environment
rm(list=ls())
# Print version of R
cat("This script was run with:", version[['version.string']], "\n")
## This script was run with: R version 4.1.2 (2021-11-01)
# Set working directory to location where script is stored
setwd(dirname(rstudioapi::getActiveDocumentContext()$path)) # requires installation of package "rstudioapi"
# Loading required libraries
require(BAS)
cat("BAS version", getNamespaceVersion("BAS"), "\n")
## BAS version 1.6.4
require(boral)
cat("boral version", getNamespaceVersion("boral"), "\n")
## boral version 2.0
require(car)
cat("car version", getNamespaceVersion("car"), "\n")
## car version 3.1-2
require(corrplot)
cat("corrplot version", getNamespaceVersion("corrplot"), "\n")
## corrplot version 0.92
require(ggplot2)
cat("ggplot2 version", getNamespaceVersion("ggplot2"), "\n")
## ggplot2 version 3.4.2
require(gplots)
cat("gplots version", getNamespaceVersion("gplots"), "\n")
## gplots version 3.1.3
require(vegan)
cat("vegan version", getNamespaceVersion("vegan"), "\n")
## vegan version 2.6-4
require(factoextra)
cat("factoextra version", getNamespaceVersion("factoextra"), "\n")
## factoextra version 1.0.7.999
require(piecewiseSEM)
cat("piecewiseSEM version", getNamespaceVersion("piecewiseSEM"), "\n")
## piecewiseSEM version 2.1.0
require(remotes)
cat("remotes version", getNamespaceVersion("remotes"), "\n")
## remotes version 2.4.2
require(nlme)
cat("nlme version", getNamespaceVersion("nlme"), "\n")
## nlme version 3.1-155
Fish were euthanized with a lethal dose of MS222 on the day of capture, following directions of the KU Leuven Animal Ethics Commission, and stored at -20 °C. Fish were kept in separate containers per site at all times. Lab based parasite screening of thawed fish involved placing individual fish in 5 or 10 ml cryo-tubes with 1 or 2 ml of distilled water. Following a vigorous shake of 10 s, the liquid was poured into a Petri dish and ectoparasites were identified and counted using a stereomicroscope. Fish were rinsed and checked again for the presence of ectoparasites on skin and fins. The intestines were examined for endoparasites. Before dissection, fish weight (± 0.1 mg) and standard length (± 1 mm) were recorded. To quantify body condition, we calculated the scaled mass index (SMI; Maceda-Veiga et al., 2014; Peig & Green, 2009). Sex was determined during dissection by inspection of gonad development. A total of 668 fish were dissected, which amounts to approximately 20 fish per location, with the exception of seven locations where only 10 fish were screened for the presence of macroparasites. Ecto- and endoparasites were morphologically identified to species level whenever possible.
# Parasite data
data <- read.csv("data_2016_2303.csv", sep=';')
data$site <- as.factor(data$site)
# Calculate parasite parameters
#names(data)
# Parasite data is overdispersed (mostly so for Trichodina), if using average abundance data, species matrix needs to be transformed
# Remove individual fish for which no parasites were counted
datao <- na.omit(data[,c(1,22:24,26:32)])
ddata <- dispweight(datao[,-1]) # Correct for overdispersion of the parasite count data
avab <- aggregate(ddata, by = list(datao[, "site"]), function(x){mean(x, na.rm =T)}) # Calculate average abundances per site
prev = aggregate(data[,c(22:24,26:32)], by = list(data[, "site"]), function(x){sum(x >0, na.rm = T)/length(x)}) # Calculate prevalence per site
medin = aggregate(data[,c(22:24,26:32)], by = list(data[, "site"]), function(x){median(x[x >0], na.rm = T)}) # Calculate median infection intensity per site
pa = aggregate(data[,c(22:24,26:32)], by = list(data[, "site"]), function(x){ifelse(mean(x, na.rm =T)>0,1,0)}) # Assess presence/absence per site
avab[is.na(avab)] <- 0
prev[is.na(prev)] <- 0
medin[is.na(medin)] <- 0
# Host condition data
avcondition <- aggregate(data$SMI, by = list(data[,1]), function(x){mean(x, na.rm =T)})[,2] # Calculate average host condition (SMI) per site
avlength <- aggregate(data$length, by = list(data[,1]), function(x){mean(x, na.rm =T)})[,2] # Calculate average host length per site
# Parasite index
sgyr <- 1:nrow(data)
stri <- 1:nrow(data)
sglu <- 1:nrow(data)
scon <- 1:nrow(data)
scysl <- 1:nrow(data)
spro <- 1:nrow(data)
saca <- 1:nrow(data)
scam <- 1:nrow(data)
sang <- 1:nrow(data)
for(j in 1:nrow(data)){
sgyr[j] <- data$Gyr[j]/sd(data$Gyr, na.rm=T)
stri[j] <- data$Tri[j]/sd(data$Tri, na.rm=T)
sglu[j] <- data$Glu[j]/sd(data$Glu, na.rm=T)
scon[j] <- data$Con[j]/sd(data$Con, na.rm=T)
scysl[j] <- data$CysL[j]/sd(data$CysL, na.rm=T)
spro[j] <- data$Pro[j]/sd(data$Pro, na.rm=T)
saca[j] <- data$Aca[j]/sd(data$Aca, na.rm=T)
scam[j] <- data$Cam[j]/sd(data$Cam, na.rm=T)
sang[j] <- data$Ang[j]/sd(data$Ang, na.rm=T)
}
PI <- 1:nrow(data)
for(j in 1:nrow(data)){
PI[j] <- 10/max(sgyr, na.rm=T)*sgyr[j] + 10/max(stri, na.rm=T)*stri[j] + 10/max(sglu, na.rm=T)*sglu[j] +
10/max(scon, na.rm=T)*scon[j] + 10/max(scysl, na.rm=T)*scysl[j] + 10/max(spro, na.rm=T)*spro[j] +
10/max(saca, na.rm=T)*saca[j] + 10/max(scam, na.rm=T)*scam[j] + 10/max(sang, na.rm=T)*sang[j]
}
PI_ecto <- 1:nrow(data)
for(j in 1:nrow(data)){
PI_ecto[j] <- 10/max(sgyr, na.rm=T)*sgyr[j] + 10/max(stri, na.rm=T)*stri[j] + 10/max(sglu, na.rm=T)*sglu[j]
}
PI_endo <- 1:nrow(data)
for(j in 1:nrow(data)){
PI_endo[j] <- 10/max(scon, na.rm=T)*scon[j] + 10/max(scysl, na.rm=T)*scysl[j] + 10/max(spro, na.rm=T)*spro[j] +
10/max(saca, na.rm=T)*saca[j] + 10/max(scam, na.rm=T)*scam[j] + 10/max(sang, na.rm=T)*sang[j]
}
avPI <- aggregate(PI, by = list(data[,1]), function(x){mean(x, na.rm =T)})[,2] # Calculate average parasite index per site
avPI_ecto <- aggregate(PI_ecto, by = list(data[,1]), function(x){mean(x, na.rm =T)})[,2] # Calculate average parasite index for ectoparasites per site
avPI_endo <- aggregate(PI_endo, by = list(data[,1]), function(x){mean(x, na.rm =T)})[,2] # Calculate average parasite index for endoparasites per site
write.csv(prev, "supplementary_table_2.csv", row.names=FALSE)
Physico-chemical data was provided by the Flemish Environmental Agency (VMM). Each fish sampling site was chosen at or near an environmental monitoring site of VMM. Water parameters include water temperature, pH, conductivity, dissolved oxygen (O2), saturation with dissolved oxygen, and Biochemical and Chemical Oxygen Demand (BOD and COD). Nutrient related water parameters include measurements of nitrate (NO3-), nitrite (NO2-), Kjeldahl nitrogen (KjN), total nitrogen (Nt), ammonium (NH4+), and total phosphorus (Pt). Following removal of strong collinear variables (significant correlation with P < 0.05 and Pearson correlation coefficient > |0.6|; Dormann et al., 2013), six environmental physico-chemical variables were retained (temperature, conductivity, COD, saturation with dissolved oxygen, ammonium, and total nitrogen), representing different aspects of water quality. For each water parameter, the average value of the year before sampling was calculated based on monthly monitoring data. Additionally, two hydromorphological variables were included: Tthe presence of a pool-riffle pattern and meanders were noted during field sampling and these parameters were included as binary variables (presence/absence) for a representation of abiotic habitat structure. Spatial (waterway) distances were calculated using the Network Analyst toolbox in ArcGIS. Upstream distance was defined as the maximal upstream distance from the sampling location. Network peripherality was calculated as the average waterway distance of a sampling location to all other locations. Hence, a total of eight environmental and two spatial variables were included in the statistical analysis.
# Environmental data (VMM)
environment <- read.csv("Environment_update.csv", sep=';')
# Spatial variables: network centrality and upstream distance
spavar <- read.csv("space2.csv", sep=';')
#plot(spavar$netcen); plot(density(spavar$netcen))
#plot(spavar$updist); plot(density(spavar$updist))
# Environmental data (from field observations)
field_data <- read.csv("field_data.csv", sep=',')
environment2 <- cbind(environment[,c(1,49,52:53,55,57,60,63)], field_data[-c(8,10,25,27,37),33:34], spavar[,c(2,3)])
environment2$pool_riffle <- as.factor(environment2$pool_riffle)
environment2$meander <- as.factor(environment2$meander)
netcen <- spavar$netcen
updist <- spavar$updist
supplementary_table_1 = cbind(environment2, netcen, updist)
write.csv(supplementary_table_1, "supplementary_table_1.csv", row.names=FALSE)
We used univariate generalized linear models to investigate how landscape-level effects modify infection patterns of individual parasite taxa, host size and condition. We kept the statistical models linear (as opposed to polynomial) and only considered main effects (i.e. no interaction terms) because we had no prior information from this study system that more complex models were necessary and because the study design with (only) 37 sampling sites was not intended for non-linear models. Ten explanatory variables (temperature, conductivity, COD, saturation with dissolved oxygen, ammonium, total nitrogen, the presence of pool-riffle patterns and meanders, upstream distance, and network peripherality) were included.
Univariate analyses - We used generalized linear models in a BMA approach to understand how infection with individual parasite taxa relate to host characteristics (length and condition), environmental conditions as well as spatial distance among sampling sites. Parasite infection was calculated in three ways at the host population level: average abundance (mean parasites per host), prevalence (percentage of infected hosts) and median infection intensity (median number of parasites in infected hosts). We calculated the individual parasitation index (IPI) following Kalbe et al. (2002) as a measurement for total parasite abundance and species richness for each individual fish. This index was calculated for all parasite species combined, and for ecto- and endoparasite species separately. For these models, we assumed a normal error distribution (which appeared to be justified, see Supplementary Figures S1-S2) and applied a Jeffrey-Zellner-Siow prior. Model assumptions (homoscedasticity of the variances and normal distribution of the errors) were assessed using the generic model plot function in R and did onlynot clearly deviate in any of the models for rare parasites. We followed a normal distribution, and not Poisson or negative binomial, for the parasite data for the common species (Trichodina sp. and Gyrodactylus spp.) and the individual parasitation index as the parameters used are deviates from count data. Rare parasites (Glugea, Contracaecum, Anguillicoloides, and unidentified cysts) were excluded from the univariate analysis because there was not enough data to obtain a good fit of the models.For rare parasites (Contracaecum sp. and Anguillicoloides crassus), we used population-level presence-absence data assuming a binomial error distribution and a uniformly distributed BIC prior. Due to low prevalences, the other parasites were not included in the species-specific models. Explanatory variables were considered important when they had a posterior inclusion probability (PIP) of 0.5. To account for overdispersion in the parasite counts, we transformed the data by downweighting overdispersed taxa following Clarke et al. (2006) using the dispweight function in the R package vegan v2.5.6 (Oksanen et al., 2013).
# Make a matrix for PIP (Posterior Inclusion Probability)
PIP <- matrix(nrow=12, ncol=14)
rownames(PIP) <- c("Host condition", "Host length", "Temperature", "Oxygen saturation", "Conductivity", "COD", "Ammonium", "Total nitrogen", "Pool riffle pattern", "Meander", "Network peripherality", "Upstream distance")
colnames(PIP) <- c("Host condition", "Host size", "Gyrodactylus abundance", "Gyrodactylus prevalence", "Gyrodactylus infection intensity", "Trichodina abundance", "Trichodina prevalence", "Trichodina infection intensity", "Glugea", "Contracaecum", "Aguillicola",
"PI", "PI ecto", "PI endo")
#Condition
bas.model <- bas.lm(avcondition ~ T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, prior="JZS")
plot(bas.model)
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(3:12),1] <- pip[2:11,1]*sign(coef.model$postmean[2:11])
#Length
bas.model <- bas.lm(avlength ~ T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(3:12),2] <- pip[2:11,1]*sign(coef.model$postmean[2:11])
#Gyrodactylus abundance
bas.model <- bas.lm(avab$Gyr ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),3] <- pip[2:13,1]*sign(coef.model$postmean[2:13])
#Gyrodactylus prevalence
bas.model <- bas.lm(prev$Gyr ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),4] <- pip[2:13,1]*sign(coef.model$postmean[2:13])
#Gyrodactylus infection intensity
bas.model <- bas.lm(medin$Gyr ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),5] <- pip[2:13,1]*sign(coef.model$postmean[2:13])
#Trichodina abundance
bas.model <- bas.lm(avab$Tri ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),6] <- pip[2:13,1]*sign(coef.model$postmean[2:13])
#Trichodina prevalence
bas.model <- bas.lm(prev$Tri ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),7] <- pip[2:13,1]*sign(coef.model$postmean[2:13])
#Trichodina infection intensity
bas.model <- bas.lm(medin$Tri ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),8] <- pip[2:13,1]*sign(coef.model$postmean[2:13])
#Glugea
bas.model <- bas.glm(pa$Glu ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, betaprior=g.prior(100), family=binomial)
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),9] <- pip[2:13,1]*sign(coef.model$postmean[2:13])
#Contracaecum
bas.model <- bas.glm(pa$Con ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, betaprior=g.prior(100), family=binomial)
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),10] <- pip[2:13,1]*sign(coef.model$postmean[2:13])
#Anguillicola
bas.model <- bas.glm(pa$Ang ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, betaprior=g.prior(100), family=binomial)
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),11] <- pip[2:13,1]*sign(coef.model$postmean[2:13])
#PI
bas.model <- bas.lm(avPI ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),12] <- pip[2:13,1]*sign(coef.model$postmean[2:13])
#PI ecto
bas.model <- bas.lm(avPI_ecto ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),13] <- pip[2:13,1]*sign(coef.model$postmean[2:13])
#PI endo
bas.model <- bas.lm(avPI_endo ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),14] <- pip[2:13,1]*sign(coef.model$postmean[2:13])
x = round(PIP, digits=2)
x[abs(PIP)<0.5] <- ""
x[abs(PIP)>0.5] <- "+"
heatmap.2(PIP[,-c(9,10,11)],
cellnote = x[,-c(9,10,11)],
#main = "Correlation",
notecex=1,
notecol="white",
density.info="none",
trace="none",
margins =c(10,8),
col=redblue(256),
dendrogram="both",
cexRow = 0.7,
cexCol = 0.7,
key.title = "PIP",
lhei = c(1,3),
lwid = c(0.5, 0.5),
srtCol = 45,
labCol = c("Host condition", "Host length",
expression(paste(italic("Gyrodactylus"), " spp. - abundance")),
expression(paste(italic("Gyrodactylus"), " spp. - prevalence")),
expression(paste(italic("Gyrodactylus"), " spp. - infection intensity")),
expression(paste(italic("Trichodina"), " sp. - abundance")),
expression(paste(italic("Trichodina"), " sp. - prevalence")),
expression(paste(italic("Trichodina"), " sp. - infection intensity")),
expression("I"[PI]*" all parasites"),
expression("I"[PI]*" ectoparasites"),
expression("I"[PI]*" endoparasites")) ,
#Colv="NA"
)
pdf("Figure2.pdf", width = 7.29, height = 4.5)
heatmap.2(PIP[,-c(9,10,11)],
cellnote = x[,-c(9,10,11)],
#main = "Correlation",
notecex=1,
notecol="white",
density.info="none",
trace="none",
margins =c(10,8),
col=redblue(256),
dendrogram="both",
cexRow = 0.7,
cexCol = 0.7,
key.title = "PIP",
lhei = c(1,3),
lwid = c(0.5, 0.5),
srtCol = 45,
labCol = c("Host condition", "Host length",
expression(paste(italic("Gyrodactylus"), " spp. - abundance")),
expression(paste(italic("Gyrodactylus"), " spp. - prevalence")),
expression(paste(italic("Gyrodactylus"), " spp. - infection intensity")),
expression(paste(italic("Trichodina"), " sp. - abundance")),
expression(paste(italic("Trichodina"), " sp. - prevalence")),
expression(paste(italic("Trichodina"), " sp. - infection intensity")),
expression("I"[PI]*" all parasites"),
expression("I"[PI]*" ectoparasites"),
expression("I"[PI]*" endoparasites")) ,
#Colv="NA"
)
dev.off()
## png
## 2
bas.model <- bas.lm(avcondition ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av +
pool_riffle + meander + netcen + updist,
data=environment2, prior="JZS")
plot(bas.model)
summary(bas.model)
## P(B != 0 | Y) model 1 model 2 model 3 model 4 model 5
## Intercept 1.00000000 1.0000 1.0000000 1.0000000 1.0000000 1.0000000
## T_av 0.12386945 0.0000 1.0000000 0.0000000 0.0000000 0.0000000
## O2_sat_av 0.05073157 0.0000 0.0000000 0.0000000 0.0000000 0.0000000
## Con_av 0.03919858 0.0000 0.0000000 0.0000000 0.0000000 0.0000000
## COD_av 0.05197281 0.0000 0.0000000 0.0000000 0.0000000 1.0000000
## NH4._av 0.04846220 0.0000 0.0000000 0.0000000 0.0000000 0.0000000
## Nt_av 0.05248386 0.0000 0.0000000 0.0000000 0.0000000 0.0000000
## pool_riffle1 0.06729856 0.0000 0.0000000 0.0000000 1.0000000 0.0000000
## meander1 0.06527705 0.0000 0.0000000 0.0000000 0.0000000 0.0000000
## netcen 0.06898689 0.0000 0.0000000 1.0000000 0.0000000 0.0000000
## updist 0.04496936 0.0000 0.0000000 0.0000000 0.0000000 0.0000000
## BF NA 1.0000 0.6511099 0.3462387 0.2957714 0.2212089
## PostProbs NA 0.6912 0.0450000 0.0239000 0.0204000 0.0153000
## R2 NA 0.0000 0.0906000 0.0565000 0.0478000 0.0315000
## dim NA 1.0000 2.0000000 2.0000000 2.0000000 2.0000000
## logmarg NA 0.0000 -0.4290769 -1.0606268 -1.2181684 -1.5086476
image(bas.model, rotate=F)
coef.model <- coef(bas.model)
#abs(coef.model$postmean)-2*coef.model$postsd > 0
plot(confint(coef.model, parm = 2:11))
## Warning in arrows(x[not.deg], ci[not.deg, 1], x[not.deg], ci[not.deg, 2], :
## zero-length arrow is of indeterminate angle and so skipped
## NULL
confint <- confint(coef.model, parm = 2:11)
write.table(confint, 'condition.txt', sep="\t")
pip <- summary(bas.model)
PIP[c(3:12),1] <- pip[2:11,1]*sign(coef.model$postmean[2:11])
#coef.model$postmean[2:11]
bas.model <- bas.lm(avlength ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av +
pool_riffle + meander + netcen + updist,
data=environment2, prior="JZS")
plot(bas.model)
summary(bas.model)
## P(B != 0 | Y) model 1 model 2 model 3 model 4 model 5
## Intercept 1.0000000 1.0000000 1.00000000 1.0000000 1.000000 1.0000000
## T_av 0.1914598 0.0000000 0.00000000 0.0000000 0.000000 0.0000000
## O2_sat_av 0.1407247 0.0000000 0.00000000 0.0000000 0.000000 0.0000000
## Con_av 0.4373524 0.0000000 0.00000000 1.0000000 1.000000 0.0000000
## COD_av 0.1445232 0.0000000 0.00000000 0.0000000 0.000000 0.0000000
## NH4._av 0.2026061 0.0000000 0.00000000 0.0000000 0.000000 0.0000000
## Nt_av 0.1889236 0.0000000 0.00000000 0.0000000 0.000000 0.0000000
## pool_riffle1 0.2259504 0.0000000 0.00000000 0.0000000 0.000000 1.0000000
## meander1 0.3217534 0.0000000 0.00000000 0.0000000 1.000000 0.0000000
## netcen 0.6121912 1.0000000 0.00000000 0.0000000 0.000000 1.0000000
## updist 0.1480063 0.0000000 0.00000000 0.0000000 0.000000 0.0000000
## BF NA 0.8278694 0.07277547 0.2894406 1.000000 0.6785609
## PostProbs NA 0.1582000 0.13910000 0.0553000 0.042500 0.0288000
## R2 NA 0.2306000 0.00000000 0.1819000 0.314300 0.2982000
## dim NA 2.0000000 1.00000000 2.0000000 3.000000 3.0000000
## logmarg NA 2.4314765 0.00000000 1.3805712 2.620376 2.2325953
image(bas.model, rotate=F)
coef.model <- coef(bas.model)
#abs(coef.model$postmean)-2*coef.model$postsd > 0
plot(confint(coef.model, parm = 2:11))
## Warning in arrows(x[not.deg], ci[not.deg, 1], x[not.deg], ci[not.deg, 2], :
## zero-length arrow is of indeterminate angle and so skipped
## Warning in arrows(x[not.deg], ci[not.deg, 1], x[not.deg], ci[not.deg, 2], :
## zero-length arrow is of indeterminate angle and so skipped
## NULL
confint <- confint(coef.model, parm = 2:11)
write.table(confint, 'length.txt', sep="\t")
pip <- summary(bas.model)
PIP[c(3:12),1] <- pip[2:11,1]*sign(coef.model$postmean[2:11])
#coef.model$postmean[2:11]
# Prediction plot
newdata = as.data.frame(cbind(rep(mean(environment2$T_av), 37),
rep(mean(environment2$O2_sat_av), 37),
rep(mean(environment2$Con_av), 37),
rep(mean(environment2$COD_av), 37),
rep(mean(environment2$NH4._av), 37),
rep(mean(environment2$Nt_av), 37),
rep(1, 37),
rep(1, 37),
rep(mean(netcen), 37),
rep(mean(updist), 37)))
colnames(newdata) <- c("T_av", "O2_sat_av", "Con_av", "COD_av", "NH4._av", "Nt_av", "pool_riffle", "meander", "netcen", "updist")
newdata[,"pool_riffle"] <- as.factor(newdata[,"pool_riffle"]); newdata[,"meander"] <- as.factor(newdata[,"meander"])
newdata1 <- newdata; newdata1[,"netcen"] <- netcen
BMA_avlength_netcen <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3j <- ggplot(environment2, aes(netcen, BMA_avlength_netcen$fit)) +
theme_bw() +
geom_line(color="red", size=1) +
geom_ribbon(aes(ymin = (BMA_avlength_netcen$fit-BMA_avlength_netcen$se.bma.fit), ymax = (BMA_avlength_netcen$fit+BMA_avlength_netcen$se.bma.fit)), alpha = .1) +
geom_point(data = environment2, aes(x=netcen, y=avlength)) +
labs(x=expression("Network peripherality [m]"), y=expression("Average host length [mm]")) +
theme(axis.title.x = element_text(size=12),
axis.title.y = element_text(size=12)) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
labs(subtitle = "j)")
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
figure3j
bas.model <- bas.lm(avab$Gyr ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, prior="JZS")
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)
summary(bas.model)
## P(B != 0 | Y) model 1 model 2 model 3 model 4 model 5
## Intercept 1.0000000 1.000000 1.00000000 1.0000000 1.0000000 1.0000000
## avlength 0.5500651 1.000000 0.00000000 0.0000000 0.0000000 1.0000000
## avcondition 0.3049274 0.000000 0.00000000 0.0000000 0.0000000 1.0000000
## T_av 0.1620852 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## O2_sat_av 0.2017699 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## Con_av 0.2621665 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## COD_av 0.8086300 1.000000 0.00000000 1.0000000 1.0000000 1.0000000
## NH4._av 0.2417809 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## Nt_av 0.8618345 1.000000 1.00000000 1.0000000 1.0000000 1.0000000
## pool_riffle1 0.1760484 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## meander1 0.1775842 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## netcen 0.3006370 0.000000 0.00000000 1.0000000 0.0000000 0.0000000
## updist 0.1648814 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## BF NA 1.000000 0.03902505 0.5452757 0.1537678 0.8271829
## PostProbs NA 0.072100 0.05160000 0.0393000 0.0369000 0.0265000
## R2 NA 0.524000 0.28790000 0.5059000 0.4100000 0.5631000
## dim NA 4.000000 2.00000000 4.0000000 3.0000000 5.0000000
## logmarg NA 6.997337 3.75378552 6.3908733 5.1250253 6.8076077
image(bas.model, rotate=F)
coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0
## [1] TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE
confint(coef.model)
## 2.5% 97.5% beta
## Intercept 9.498176e-02 2.044717e-01 1.489134e-01
## avlength -2.515130e-02 2.264807e-05 -8.122675e-03
## avcondition -9.805559e-01 0.000000e+00 -1.500326e-01
## T_av -2.205915e-02 1.634873e-02 -1.493132e-04
## O2_sat_av -1.678010e-03 5.033937e-03 3.794809e-04
## Con_av -9.556692e-05 5.266581e-04 5.726233e-05
## COD_av 0.000000e+00 1.300051e-02 6.652983e-03
## NH4._av -6.709599e-02 7.179572e-02 -1.888173e-04
## Nt_av 0.000000e+00 5.230759e-02 2.930918e-02
## pool_riffle1 -3.796033e-02 1.029158e-01 3.966550e-03
## meander1 -9.234735e-02 5.004225e-02 -4.806259e-03
## netcen -5.720044e-07 1.018422e-05 1.332021e-06
## updist -1.325467e-06 1.732882e-06 4.618571e-08
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"
plot(confint(coef.model, parm = 2:11))
## NULL
confint <- confint(coef.model, parm = 2:11)
write.table(confint, 'GyroAA.txt', sep="\t")
pip <- summary(bas.model)
PIP[c(1:12),3] <- pip[2:13,1]*sign(coef.model$postmean[2:13])
# Prediction plot
newdata = as.data.frame(cbind(rep(mean(avlength), 37),
rep(mean(avcondition), 37),
rep(mean(environment2$T_av), 37),
rep(mean(environment2$O2_sat_av), 37),
rep(mean(environment2$Con_av), 37),
rep(mean(environment2$COD_av), 37),
rep(mean(environment2$NH4._av), 37),
rep(mean(environment2$Nt_av), 37),
rep(1, 37),
rep(1, 37),
rep(mean(netcen), 37),
rep(mean(updist), 37)))
colnames(newdata) <- c("avlength", "avcondition", "T_av", "O2_sat_av", "Con_av", "COD_av", "NH4._av", "Nt_av", "pool_riffle", "meander", "netcen", "updist")
newdata[,"pool_riffle"] <- as.factor(newdata[,"pool_riffle"]); newdata[,"meander"] <- as.factor(newdata[,"meander"])
newdata1 <- newdata; newdata1[,"avcondition"] <- avcondition
BMA_Gyr_avcond <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3k = ggplot(environment2, aes(avcondition, BMA_Gyr_avcond$fit)) +
theme_bw() +
geom_line(color="red", size=1) +
geom_ribbon(aes(ymin = (BMA_Gyr_avcond$fit-BMA_Gyr_avcond$se.bma.fit), ymax = (BMA_Gyr_avcond$fit+BMA_Gyr_avcond$se.bma.fit)), alpha = .1) +
geom_point(data = environment2, aes(x=avcondition, y=avab$Gyr)) +
labs(x=expression("Average host condition"), y=expression(paste(italic("Gyrodactylus"), " spp. - abundance"))) +
theme(axis.title.x = element_text(size=12),
axis.title.y = element_text(size=12)) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
labs(subtitle = "k)")
figure3k
newdata1 <- newdata; newdata1[,"COD_av"] <- environment2$COD_av
BMA_Gyr_COD_av <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3a = ggplot(environment2, aes(COD_av, BMA_Gyr_COD_av$fit)) +
theme_bw() +
geom_line(color="red", size=1) +
geom_ribbon(aes(ymin = (BMA_Gyr_COD_av$fit-BMA_Gyr_COD_av$se.bma.fit), ymax = (BMA_Gyr_COD_av$fit+BMA_Gyr_COD_av$se.bma.fit)), alpha = .1) +
geom_point(data = environment2, aes(x=COD_av, y=avab$Gyr)) +
labs(x=expression("Chemical Oxygen Demand [mg/L]"), y=expression(paste(italic("Gyrodactylus"), " spp. - abundance"))) +
theme(axis.title.x = element_text(size=12),
axis.title.y = element_text(size=12)) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
labs(subtitle = "a)")
figure3a
newdata1 <- newdata; newdata1[,"Nt_av"] <- environment2$Nt_av
BMA_Gyr_Nt_av <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3e = ggplot(environment2, aes(Nt_av, BMA_Gyr_Nt_av$fit)) +
theme_bw() +
geom_line(color="red", size=1) +
geom_ribbon(aes(ymin = (BMA_Gyr_Nt_av$fit-BMA_Gyr_Nt_av$se.bma.fit), ymax = (BMA_Gyr_Nt_av$fit+BMA_Gyr_Nt_av$se.bma.fit)), alpha = .1) +
geom_point(data = environment2, aes(x=Nt_av, y=avab$Gyr)) +
labs(x=expression("Nt [mg/L]"), y=expression(paste(italic("Gyrodactylus"), " spp. - abundance"))) +
theme(axis.title.x = element_text(size=12),
axis.title.y = element_text(size=12)) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
labs(subtitle = "e)")
figure3e
bas.model <- bas.lm(medin$Gyr ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, prior="JZS")
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)
summary(bas.model)
## P(B != 0 | Y) model 1 model 2 model 3 model 4 model 5
## Intercept 1.00000000 1.0000 1.0000000 1.000000 1.0000000 1.0000000
## avlength 0.02065469 0.0000 0.0000000 0.000000 0.0000000 0.0000000
## avcondition 0.04286905 0.0000 1.0000000 0.000000 0.0000000 0.0000000
## T_av 0.02991678 0.0000 0.0000000 0.000000 0.0000000 0.0000000
## O2_sat_av 0.02400544 0.0000 0.0000000 0.000000 0.0000000 0.0000000
## Con_av 0.02143314 0.0000 0.0000000 0.000000 0.0000000 0.0000000
## COD_av 0.03655302 0.0000 0.0000000 0.000000 1.0000000 0.0000000
## NH4._av 0.03250611 0.0000 0.0000000 0.000000 0.0000000 1.0000000
## Nt_av 0.04088810 0.0000 0.0000000 1.000000 0.0000000 0.0000000
## pool_riffle1 0.02307275 0.0000 0.0000000 0.000000 0.0000000 0.0000000
## meander1 0.02038443 0.0000 0.0000000 0.000000 0.0000000 0.0000000
## netcen 0.02173649 0.0000 0.0000000 0.000000 0.0000000 0.0000000
## updist 0.02738794 0.0000 0.0000000 0.000000 0.0000000 0.0000000
## BF NA 1.0000 0.3176354 0.302852 0.2488803 0.2411373
## PostProbs NA 0.7775 0.0206000 0.019600 0.0161000 0.0156000
## R2 NA 0.0000 0.0517000 0.049100 0.0381000 0.0363000
## dim NA 1.0000 2.0000000 2.000000 2.0000000 2.0000000
## logmarg NA 0.0000 -1.1468512 -1.194511 -1.3907833 -1.4223890
image(bas.model, rotate=F)
coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0
## [1] TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE
confint(coef.model)
## 2.5% 97.5% beta
## Intercept 1.564588 3.19689 2.405405e+00
## avlength 0.000000 0.00000 -5.938396e-04
## avcondition 0.000000 0.00000 -2.939823e-01
## T_av 0.000000 0.00000 -7.368692e-03
## O2_sat_av 0.000000 0.00000 -4.581964e-04
## Con_av 0.000000 0.00000 -1.594570e-05
## COD_av 0.000000 0.00000 1.609035e-03
## NH4._av 0.000000 0.00000 1.027436e-02
## Nt_av 0.000000 0.00000 7.570653e-03
## pool_riffle1 0.000000 0.00000 1.051497e-02
## meander1 0.000000 0.00000 2.511911e-03
## netcen 0.000000 0.00000 3.990331e-07
## updist 0.000000 0.00000 4.560731e-07
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"
plot(confint(coef.model, parm = 2:11))
## NULL
confint <- confint(coef.model, parm = 2:11)
write.table(confint, 'GyroAA.txt', sep="\t")
pip <- summary(bas.model)
PIP[c(1:12),3] <- pip[2:13,1]*sign(coef.model$postmean[2:13])
bas.model <- bas.lm(prev$Gyr ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, prior="JZS")
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)
summary(bas.model)
## P(B != 0 | Y) model 1 model 2 model 3 model 4 model 5
## Intercept 1.00000000 1.000000 1.00000000 1.0000000 1.0000000 1.00000000
## avlength 0.16685753 0.000000 0.00000000 0.0000000 1.0000000 0.00000000
## avcondition 0.07728991 0.000000 0.00000000 0.0000000 0.0000000 0.00000000
## T_av 0.07242807 0.000000 0.00000000 0.0000000 0.0000000 0.00000000
## O2_sat_av 0.09667222 0.000000 0.00000000 0.0000000 0.0000000 0.00000000
## Con_av 0.58076129 1.000000 0.00000000 1.0000000 0.0000000 0.00000000
## COD_av 0.17268847 0.000000 0.00000000 1.0000000 0.0000000 0.00000000
## NH4._av 0.09134534 0.000000 0.00000000 0.0000000 0.0000000 0.00000000
## Nt_av 0.11056039 0.000000 0.00000000 0.0000000 0.0000000 0.00000000
## pool_riffle1 0.13922799 0.000000 0.00000000 0.0000000 0.0000000 0.00000000
## meander1 0.09905733 0.000000 0.00000000 0.0000000 0.0000000 0.00000000
## netcen 0.11864979 0.000000 0.00000000 0.0000000 0.0000000 1.00000000
## updist 0.06837623 0.000000 0.00000000 0.0000000 0.0000000 0.00000000
## BF NA 1.000000 0.08283414 0.8857497 0.1246735 0.07409002
## PostProbs NA 0.228800 0.22750000 0.0369000 0.0285000 0.01700000
## R2 NA 0.233300 0.00000000 0.3039000 0.1341000 0.10730000
## dim NA 2.000000 1.00000000 3.0000000 2.0000000 2.00000000
## logmarg NA 2.490915 0.00000000 2.3695941 0.4088580 -0.11155941
image(bas.model, rotate=F)
coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0
## [1] TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE
confint(coef.model)
## 2.5% 97.5% beta
## Intercept 2.498255e-01 3.931726e-01 3.224082e-01
## avlength -2.252212e-02 2.446193e-04 -2.431148e-03
## avcondition -4.252456e-01 6.599841e-03 -1.986363e-02
## T_av -1.220894e-02 6.769487e-03 -7.123130e-04
## O2_sat_av -4.258310e-03 6.559308e-05 -2.662986e-04
## Con_av 0.000000e+00 8.103246e-04 3.112530e-04
## COD_av 0.000000e+00 7.914329e-03 9.376253e-04
## NH4._av -1.771781e-03 4.657468e-02 1.064409e-03
## Nt_av 0.000000e+00 2.507558e-02 1.612509e-03
## pool_riffle1 0.000000e+00 1.573096e-01 1.465705e-02
## meander1 -1.017710e-01 1.415960e-02 -6.966751e-03
## netcen -1.412903e-07 7.995999e-06 5.765326e-07
## updist -7.999739e-07 6.404386e-07 -6.922933e-09
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"
plot(confint(coef.model, parm = 2:11))
## NULL
confint <- confint(coef.model, parm = 2:11)
write.table(confint, 'GyroAA.txt', sep="\t")
pip <- summary(bas.model)
PIP[c(1:12),3] <- pip[2:13,1]*sign(coef.model$postmean[2:13])
newdata1 <- newdata; newdata1[,"Con_av"] <- environment2$Con_av
BMA_Gyr_Con_av <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3i = ggplot(environment2, aes(Con_av, BMA_Gyr_Con_av$fit)) +
theme_bw() +
geom_line(color="red", size=1) +
geom_ribbon(aes(ymin = (BMA_Gyr_Con_av$fit-BMA_Gyr_Con_av$se.bma.fit), ymax = (BMA_Gyr_Con_av$fit+BMA_Gyr_Con_av$se.bma.fit)), alpha = .1) +
geom_point(data = environment2, aes(x=Con_av, y=prev$Gyr)) +
labs(x=expression(paste("Conductivity [", mu, "S/cm]")), y=expression(paste(italic("Gyrodactylus"), " spp. - prevalence"))) +
theme(axis.title.x = element_text(size=12),
axis.title.y = element_text(size=12)) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
labs(subtitle = "i)")
figure3i
bas.model <- bas.lm(avab$Tri ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, prior="JZS")
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)
summary(bas.model)
## P(B != 0 | Y) model 1 model 2 model 3 model 4 model 5
## Intercept 1.00000000 1.0000000 1.000000 1.000000 1.0000000 1.0000000
## avlength 0.12444042 0.0000000 0.000000 0.000000 0.0000000 0.0000000
## avcondition 0.09559152 0.0000000 0.000000 0.000000 0.0000000 0.0000000
## T_av 0.08500185 0.0000000 0.000000 0.000000 0.0000000 0.0000000
## O2_sat_av 0.09751169 0.0000000 0.000000 0.000000 0.0000000 0.0000000
## Con_av 0.42180029 0.0000000 1.000000 0.000000 0.0000000 0.0000000
## COD_av 0.46724539 0.0000000 1.000000 0.000000 0.0000000 1.0000000
## NH4._av 0.21465122 0.0000000 0.000000 0.000000 1.0000000 0.0000000
## Nt_av 0.34576751 0.0000000 0.000000 1.000000 0.0000000 1.0000000
## pool_riffle1 0.13743699 0.0000000 0.000000 0.000000 0.0000000 0.0000000
## meander1 0.08968689 0.0000000 0.000000 0.000000 0.0000000 0.0000000
## netcen 0.10313787 0.0000000 0.000000 0.000000 0.0000000 0.0000000
## updist 0.08604777 0.0000000 0.000000 0.000000 0.0000000 0.0000000
## BF NA 0.0239404 1.000000 0.170643 0.1105942 0.2711811
## PostProbs NA 0.1650000 0.104400 0.098000 0.0635000 0.0283000
## R2 NA 0.0000000 0.358500 0.209300 0.1890000 0.3063000
## dim NA 1.0000000 3.000000 2.000000 2.0000000 3.0000000
## logmarg NA 0.0000000 3.732188 1.964006 1.5303004 2.4272192
image(bas.model, rotate=F)
coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0
## [1] TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE
confint(coef.model)
## 2.5% 97.5% beta
## Intercept 1.266097e-01 2.732195e-01 2.002413e-01
## avlength -1.608394e-02 3.767937e-04 -1.114390e-03
## avcondition -4.632070e-01 1.274208e-01 -2.394139e-02
## T_av -1.438322e-02 1.013193e-02 -1.479738e-04
## O2_sat_av -1.905880e-03 3.300222e-03 1.132458e-04
## Con_av 0.000000e+00 7.774852e-04 2.109076e-04
## COD_av 0.000000e+00 1.393205e-02 4.203136e-03
## NH4._av -2.205558e-03 1.033887e-01 1.096330e-02
## Nt_av -4.771345e-05 4.961465e-02 1.044647e-02
## pool_riffle1 0.000000e+00 1.466687e-01 1.212839e-02
## meander1 -9.597451e-02 4.095486e-03 -1.678368e-03
## netcen -6.164040e-06 1.239477e-06 -2.846749e-07
## updist -1.511702e-06 6.054639e-07 -1.915234e-08
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"
plot(confint(coef.model, parm = 2:11))
## NULL
confint <- confint(coef.model, parm = 2:11)
write.table(confint, 'GyroAA.txt', sep="\t")
pip <- summary(bas.model)
PIP[c(1:12),3] <- pip[2:13,1]*sign(coef.model$postmean[2:13])
bas.model <- bas.lm(medin$Tri ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, prior="JZS")
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)
summary(bas.model)
## P(B != 0 | Y) model 1 model 2 model 3 model 4 model 5
## Intercept 1.0000000 1.000000 1.00000000 1.0000000 1.0000000 1.0000000
## avlength 0.2881774 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## avcondition 0.2224163 0.000000 0.00000000 0.0000000 1.0000000 0.0000000
## T_av 0.1574615 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## O2_sat_av 0.1621652 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## Con_av 0.7912882 1.000000 0.00000000 1.0000000 1.0000000 1.0000000
## COD_av 0.7219333 1.000000 0.00000000 0.0000000 1.0000000 1.0000000
## NH4._av 0.3248706 0.000000 1.00000000 1.0000000 0.0000000 0.0000000
## Nt_av 0.3218716 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## pool_riffle1 0.2081402 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## meander1 0.2655250 0.000000 0.00000000 0.0000000 0.0000000 1.0000000
## netcen 0.2107764 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## updist 0.1536937 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## BF NA 1.000000 0.04957628 0.2279222 0.5092549 0.4837322
## PostProbs NA 0.148800 0.04060000 0.0339000 0.0227000 0.0216000
## R2 NA 0.449900 0.26810000 0.3987000 0.4816000 0.4799000
## dim NA 3.000000 2.00000000 3.0000000 4.0000000 4.0000000
## logmarg NA 6.288599 3.28435583 4.8098478 5.6137920 5.5623749
image(bas.model, rotate=F)
coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0
## [1] TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE
confint(coef.model)
## 2.5% 97.5% beta
## Intercept 1.660851e+01 3.395723e+01 2.547297e+01
## avlength -3.541722e+00 2.164229e-02 -5.158998e-01
## avcondition -1.441975e+02 4.413419e+00 -1.375058e+01
## T_av -1.880970e+00 4.198759e+00 2.065105e-01
## O2_sat_av -5.817500e-01 4.198218e-01 -1.498026e-02
## Con_av 0.000000e+00 1.277170e-01 6.585727e-02
## COD_av 0.000000e+00 1.980322e+00 9.195418e-01
## NH4._av -1.801730e+00 1.628155e+01 2.222211e+00
## Nt_av -1.003826e-01 6.309924e+00 1.012166e+00
## pool_riffle1 -4.799011e-01 2.544697e+01 2.247962e+00
## meander1 -2.640241e+01 1.529459e-01 -3.523477e+00
## netcen -1.457582e-03 1.774935e-04 -1.366046e-04
## updist -2.281233e-04 2.696601e-04 5.653370e-06
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"
plot(confint(coef.model, parm = 2:11))
## NULL
confint <- confint(coef.model, parm = 2:11)
write.table(confint, 'GyroAA.txt', sep="\t")
pip <- summary(bas.model)
PIP[c(1:12),3] <- pip[2:13,1]*sign(coef.model$postmean[2:13])
newdata1 <- newdata; newdata1[,"Con_av"] <- environment2$Con_av
BMA_Tri_Con_av <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3h = ggplot(environment2, aes(Con_av, BMA_Tri_Con_av$fit)) +
theme_bw() +
geom_line(color="red", size=1) +
geom_ribbon(aes(ymin = (BMA_Tri_Con_av$fit-BMA_Tri_Con_av$se.bma.fit), ymax = (BMA_Tri_Con_av$fit+BMA_Tri_Con_av$se.bma.fit)), alpha = .1) +
geom_point(data = environment2, aes(x=Con_av, y=medin$Tri)) +
labs(x=expression(paste("Conductivity [", mu, "S/cm]")), y=expression(paste(italic("Trichodina"), " sp. - infection intensity"))) +
theme(axis.title.x = element_text(size=12),
axis.title.y = element_text(size=12)) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
labs(subtitle = "h)")
figure3h
newdata1 <- newdata; newdata1[,"COD_av"] <- environment2$COD_av
BMA_Tri_COD_av <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3d = ggplot(environment2, aes(COD_av, BMA_Tri_COD_av$fit)) +
theme_bw() +
geom_line(color="red", size=1) +
geom_ribbon(aes(ymin = (BMA_Tri_COD_av$fit-BMA_Tri_COD_av$se.bma.fit), ymax = (BMA_Tri_COD_av$fit+BMA_Tri_COD_av$se.bma.fit)), alpha = .1) +
geom_point(data = environment2, aes(x=COD_av, y=medin$Tri)) +
labs(x=expression(paste("Chemical Oxygen Demand [mg/L]")), y=expression(paste(italic("Trichodina"), " sp. - infection intensity"))) +
theme(axis.title.x = element_text(size=12),
axis.title.y = element_text(size=12)) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
labs(subtitle = "d)")
figure3d
bas.model <- bas.lm(prev$Tri ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, prior="JZS")
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)
summary(bas.model)
## P(B != 0 | Y) model 1 model 2 model 3 model 4 model 5
## Intercept 1.00000000 1.0000000 1.0000000 1.0000000 1.0000000 1.000000
## avlength 0.04005815 0.0000000 0.0000000 0.0000000 0.0000000 0.000000
## avcondition 0.03389439 0.0000000 0.0000000 0.0000000 0.0000000 0.000000
## T_av 0.03654853 0.0000000 0.0000000 0.0000000 0.0000000 0.000000
## O2_sat_av 0.03749123 0.0000000 0.0000000 0.0000000 0.0000000 0.000000
## Con_av 0.20573471 0.0000000 1.0000000 0.0000000 1.0000000 0.000000
## COD_av 0.04629842 0.0000000 0.0000000 0.0000000 0.0000000 0.000000
## NH4._av 0.06690739 0.0000000 0.0000000 1.0000000 0.0000000 0.000000
## Nt_av 0.04604656 0.0000000 0.0000000 0.0000000 0.0000000 1.000000
## pool_riffle1 0.03411173 0.0000000 0.0000000 0.0000000 0.0000000 0.000000
## meander1 0.03657952 0.0000000 0.0000000 0.0000000 0.0000000 0.000000
## netcen 0.07467299 0.0000000 0.0000000 0.0000000 1.0000000 0.000000
## updist 0.07352673 0.0000000 0.0000000 0.0000000 0.0000000 0.000000
## BF NA 0.5616705 0.7274201 0.2732671 1.0000000 0.155205
## PostProbs NA 0.6354000 0.0686000 0.0258000 0.0171000 0.014600
## R2 NA 0.0000000 0.1264000 0.0750000 0.2249000 0.044000
## dim NA 1.0000000 2.0000000 2.0000000 3.0000000 2.000000
## logmarg NA 0.0000000 0.2585887 -0.7204656 0.5768398 -1.286169
image(bas.model, rotate=F)
coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0
## [1] TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE
confint(coef.model)
## 2.5% 97.5% beta
## Intercept 7.818130e-01 8.857532e-01 8.340580e-01
## avlength 0.000000e+00 0.000000e+00 -1.453558e-04
## avcondition 0.000000e+00 0.000000e+00 -2.610917e-03
## T_av 0.000000e+00 0.000000e+00 -1.135202e-05
## O2_sat_av 0.000000e+00 0.000000e+00 -1.893466e-05
## Con_av -1.787956e-07 4.202091e-04 6.570599e-05
## COD_av 0.000000e+00 0.000000e+00 9.451636e-05
## NH4._av -3.379414e-05 1.957452e-02 1.713706e-03
## Nt_av 0.000000e+00 0.000000e+00 1.739937e-04
## pool_riffle1 0.000000e+00 0.000000e+00 2.104077e-04
## meander1 0.000000e+00 0.000000e+00 -6.652944e-04
## netcen -3.742651e-06 9.178620e-09 -3.317353e-07
## updist -1.636752e-06 5.369651e-10 -1.309710e-07
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"
plot(confint(coef.model, parm = 2:11))
## NULL
confint <- confint(coef.model, parm = 2:11)
write.table(confint, 'GyroAA.txt', sep="\t")
pip <- summary(bas.model)
PIP[c(1:12),3] <- pip[2:13,1]*sign(coef.model$postmean[2:13])
bas.model <- bas.glm(pa$Glu ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av +
NH4._av + Nt_av + SM_av + pool_riffle + meander +
spavar$netcen + spavar$updist, data=environment2, betaprior=g.prior(100), family=binomial)
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)
summary(bas.model)
## P(B != 0 | Y) model 1 model 2 model 3 model 4
## Intercept 1.0000000 1.000000 1.0000000 1.000000 1.0000000
## avlength 0.6721313 0.000000 1.0000000 0.000000 0.0000000
## avcondition 0.4472412 0.000000 1.0000000 1.000000 1.0000000
## T_av 0.2835083 0.000000 0.0000000 0.000000 0.0000000
## O2_sat_av 0.3836304 0.000000 0.0000000 0.000000 0.0000000
## Con_av 0.9992554 1.000000 1.0000000 1.000000 1.0000000
## COD_av 0.5410645 1.000000 1.0000000 0.000000 1.0000000
## NH4._av 0.8430908 1.000000 1.0000000 1.000000 1.0000000
## Nt_av 0.5500488 1.000000 0.0000000 1.000000 1.0000000
## SM_av 0.4324219 1.000000 0.0000000 1.000000 0.0000000
## pool_riffle1 0.4305176 0.000000 0.0000000 0.000000 0.0000000
## meander1 0.9996338 1.000000 1.0000000 1.000000 1.0000000
## spavar$netcen 0.4574951 0.000000 0.0000000 1.000000 0.0000000
## spavar$updist 0.8283325 1.000000 1.0000000 0.000000 1.0000000
## BF NA 1.000000 0.8189871 0.769804 0.7351301
## PostProbs NA 0.053100 0.0522000 0.039900 0.0384000
## R2 NA 1.000000 1.0000000 1.000000 1.0000000
## dim NA 8.000000 8.0000000 8.000000 8.0000000
## logmarg NA -5.808403 -6.0080902 -6.070023 -6.1161110
## model 5
## Intercept 1.0000000
## avlength 1.0000000
## avcondition 0.0000000
## T_av 0.0000000
## O2_sat_av 1.0000000
## Con_av 1.0000000
## COD_av 0.0000000
## NH4._av 0.0000000
## Nt_av 0.0000000
## SM_av 0.0000000
## pool_riffle1 1.0000000
## meander1 1.0000000
## spavar$netcen 1.0000000
## spavar$updist 1.0000000
## BF 0.5830791
## PostProbs 0.0315000
## R2 1.0000000
## dim 8.0000000
## logmarg -6.3478357
image(bas.model, rotate=F)
coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0
## [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE FALSE
confint(coef.model)
## 2.5% 97.5% beta
## Intercept -2.877566e+06 3.000550e+06 -6.056996e+03
## avlength -4.420045e+04 4.258743e+04 8.562288e+01
## avcondition -8.578579e+05 9.088818e+05 -1.218086e+01
## T_av -3.140698e+04 2.899503e+04 1.760489e+01
## O2_sat_av -5.868290e+03 5.525430e+03 8.380211e+00
## Con_av -1.518537e+03 1.619110e+03 4.195005e+00
## COD_av -1.323990e+04 1.518345e+04 9.690644e+00
## NH4._av -2.252463e+05 2.166739e+05 -2.691592e+02
## Nt_av -4.470842e+04 4.775553e+04 3.669994e+01
## SM_av -2.069908e+03 2.430891e+03 -1.161863e+00
## pool_riffle1 -1.551102e+05 1.439102e+05 8.409678e+01
## meander1 -4.719167e+05 5.411526e+05 -1.301335e+03
## spavar$netcen -8.973198e+00 1.065294e+01 -8.483558e-03
## spavar$updist -5.016680e+00 5.104948e+00 -7.734368e-03
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"
plot(confint(coef.model, parm = 2:11))
## NULL
confint <- confint(coef.model, parm = 2:11)
pip <- summary(bas.model)
PIP[c(1:12),9] <- pip[2:13,1]*sign(coef.model$postmean[2:13])
bas.model <- bas.glm(pa$Con ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av +
NH4._av + Nt_av + SM_av + pool_riffle + meander +
spavar$netcen + spavar$updist, data=environment2, betaprior=g.prior(100), family=binomial)
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)
summary(bas.model)
## P(B != 0 | Y) model 1 model 2 model 3 model 4
## Intercept 1.00000000 1.00000 1.0000000 1.0000000 1.0000000
## avlength 0.08128662 0.00000 1.0000000 0.0000000 0.0000000
## avcondition 0.02805176 0.00000 0.0000000 1.0000000 0.0000000
## T_av 0.01058350 0.00000 0.0000000 0.0000000 0.0000000
## O2_sat_av 0.02976074 0.00000 0.0000000 0.0000000 1.0000000
## Con_av 0.01105957 0.00000 0.0000000 0.0000000 0.0000000
## COD_av 0.01395264 0.00000 0.0000000 0.0000000 0.0000000
## NH4._av 0.01990967 0.00000 0.0000000 0.0000000 0.0000000
## Nt_av 0.01219482 0.00000 0.0000000 0.0000000 0.0000000
## SM_av 0.01367188 0.00000 0.0000000 0.0000000 0.0000000
## pool_riffle1 0.01322021 0.00000 0.0000000 0.0000000 0.0000000
## meander1 0.01164551 0.00000 0.0000000 0.0000000 0.0000000
## spavar$netcen 0.01945801 0.00000 0.0000000 0.0000000 0.0000000
## spavar$updist 0.01293945 0.00000 0.0000000 0.0000000 0.0000000
## BF NA 1.00000 0.9507459 0.2581488 0.2583995
## PostProbs NA 0.77910 0.0562000 0.0172000 0.0172000
## R2 NA 0.00000 0.0864000 0.0365000 0.0366000
## dim NA 1.00000 2.0000000 2.0000000 2.0000000
## logmarg NA -25.82594 -25.8764468 -27.1801577 -27.1791867
## model 5
## Intercept 1.00000
## avlength 0.00000
## avcondition 0.00000
## T_av 0.00000
## O2_sat_av 0.00000
## Con_av 0.00000
## COD_av 0.00000
## NH4._av 1.00000
## Nt_av 0.00000
## SM_av 0.00000
## pool_riffle1 0.00000
## meander1 0.00000
## spavar$netcen 0.00000
## spavar$updist 0.00000
## BF 0.19791
## PostProbs 0.01350
## R2 0.02640
## dim 2.00000
## logmarg -27.44588
image(bas.model, rotate=F)
coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0
## [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE FALSE
confint(coef.model)
## 2.5% 97.5% beta
## Intercept -9.880859 2.3742736 -5.430024e-01
## avlength 0.000000 0.1497804 1.509881e-02
## avcondition 0.000000 0.0000000 -1.814762e-01
## T_av 0.000000 0.0000000 7.945717e-04
## O2_sat_av 0.000000 0.0000000 1.286763e-03
## Con_av 0.000000 0.0000000 5.456325e-06
## COD_av 0.000000 0.0000000 -3.836691e-04
## NH4._av 0.000000 0.0000000 -6.098874e-03
## Nt_av 0.000000 0.0000000 -8.846247e-04
## SM_av 0.000000 0.0000000 1.337787e-04
## pool_riffle1 0.000000 0.0000000 -5.842378e-03
## meander1 0.000000 0.0000000 3.379776e-03
## spavar$netcen 0.000000 0.0000000 -8.679200e-07
## spavar$updist 0.000000 0.0000000 -1.196324e-07
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"
plot(confint(coef.model, parm = 2:11))
## NULL
confint <- confint(coef.model, parm = 2:11)
pip <- summary(bas.model)
PIP[c(1:12),10] <- pip[2:13,1]*sign(coef.model$postmean[2:13])
bas.model <- bas.glm(pa$Ang ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av +
NH4._av + Nt_av + SM_av + pool_riffle + meander +
spavar$netcen + spavar$updist, data=environment2, betaprior=g.prior(100), family=binomial)
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)
summary(bas.model)
## P(B != 0 | Y) model 1 model 2 model 3 model 4
## Intercept 1.0000000 1.000000 1.0000000 1.0000000 1.0000000
## avlength 0.9849854 1.000000 1.0000000 1.0000000 1.0000000
## avcondition 0.3566895 0.000000 0.0000000 0.0000000 0.0000000
## T_av 0.9871948 1.000000 1.0000000 1.0000000 1.0000000
## O2_sat_av 0.8635742 1.000000 1.0000000 1.0000000 1.0000000
## Con_av 0.5832520 0.000000 1.0000000 1.0000000 0.0000000
## COD_av 0.3989502 0.000000 0.0000000 1.0000000 0.0000000
## NH4._av 0.9910645 1.000000 1.0000000 1.0000000 1.0000000
## Nt_av 0.5583984 0.000000 1.0000000 0.0000000 0.0000000
## SM_av 0.9698608 1.000000 1.0000000 1.0000000 1.0000000
## pool_riffle1 0.4013184 0.000000 0.0000000 0.0000000 1.0000000
## meander1 0.9865601 1.000000 1.0000000 1.0000000 1.0000000
## spavar$netcen 0.8922974 1.000000 1.0000000 1.0000000 1.0000000
## spavar$updist 0.5290039 1.000000 0.0000000 0.0000000 1.0000000
## BF NA 1.000000 0.4728289 0.3794633 0.2858158
## PostProbs NA 0.080800 0.0735000 0.0537000 0.0429000
## R2 NA 1.000000 1.0000000 1.0000000 1.0000000
## dim NA 9.000000 10.0000000 10.0000000 10.0000000
## logmarg NA -9.802499 -10.5515211 -10.7714968 -11.0549070
## model 5
## Intercept 1.0000000
## avlength 1.0000000
## avcondition 1.0000000
## T_av 1.0000000
## O2_sat_av 0.0000000
## Con_av 1.0000000
## COD_av 0.0000000
## NH4._av 1.0000000
## Nt_av 1.0000000
## SM_av 1.0000000
## pool_riffle1 0.0000000
## meander1 1.0000000
## spavar$netcen 1.0000000
## spavar$updist 0.0000000
## BF 0.2037316
## PostProbs 0.0324000
## R2 1.0000000
## dim 10.0000000
## logmarg -11.3934513
image(bas.model, rotate=F)
coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0
## [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE FALSE
confint(coef.model)
## 2.5% 97.5% beta
## Intercept -2.427843e+06 2.407495e+06 -6.443834e+03
## avlength -4.661497e+04 4.211047e+04 1.321190e+02
## avcondition -8.356648e+05 8.401762e+05 -3.035671e+02
## T_av -1.531048e+05 1.561560e+05 5.865874e+02
## O2_sat_av -1.899200e+04 1.821646e+04 -3.714070e+01
## Con_av -6.375893e+02 7.129852e+02 3.717928e-01
## COD_av -7.961092e+03 7.947834e+03 1.910391e+00
## NH4._av -3.591783e+05 3.383641e+05 -1.053055e+03
## Nt_av -4.799272e+04 3.935074e+04 -4.605084e+01
## SM_av -5.743284e+03 5.697560e+03 1.280831e+01
## pool_riffle1 -1.496669e+05 1.333208e+05 -1.219107e+02
## meander1 -4.050513e+05 3.841885e+05 -1.391517e+03
## spavar$netcen -1.550874e+01 1.393704e+01 -2.782679e-02
## spavar$updist -3.150044e+00 2.889053e+00 1.171764e-03
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"
plot(confint(coef.model, parm = 2:11))
## NULL
confint <- confint(coef.model, parm = 2:11)
pip <- summary(bas.model)
PIP[c(1:12),11] <- pip[2:13,1]*sign(coef.model$postmean[2:13])
bas.model <- bas.lm(avPI ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, prior="JZS")
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)
summary(bas.model)
## P(B != 0 | Y) model 1 model 2 model 3 model 4 model 5
## Intercept 1.0000000 1.00000000 1.00000000 1.00000000 1.000000 1.000000
## avlength 0.2261960 0.00000000 0.00000000 0.00000000 0.000000 0.000000
## avcondition 0.1810890 0.00000000 0.00000000 0.00000000 0.000000 0.000000
## T_av 0.4574073 0.00000000 1.00000000 0.00000000 1.000000 0.000000
## O2_sat_av 0.1952767 0.00000000 0.00000000 0.00000000 0.000000 0.000000
## Con_av 0.4598919 0.00000000 0.00000000 0.00000000 0.000000 1.000000
## COD_av 0.5684747 0.00000000 0.00000000 0.00000000 0.000000 1.000000
## NH4._av 0.2748144 0.00000000 0.00000000 0.00000000 0.000000 0.000000
## Nt_av 0.4582802 0.00000000 0.00000000 1.00000000 1.000000 0.000000
## pool_riffle1 0.3054914 0.00000000 0.00000000 0.00000000 0.000000 0.000000
## meander1 0.1985617 0.00000000 0.00000000 0.00000000 0.000000 0.000000
## netcen 0.1910872 0.00000000 0.00000000 0.00000000 0.000000 0.000000
## updist 0.4153657 0.00000000 0.00000000 0.00000000 0.000000 1.000000
## BF NA 0.01969453 0.09251289 0.08737546 0.355023 1.000000
## PostProbs NA 0.10850000 0.04250000 0.04010000 0.029600 0.025000
## R2 NA 0.00000000 0.18980000 0.18710000 0.325400 0.424800
## dim NA 1.00000000 2.00000000 2.00000000 3.000000 4.000000
## logmarg NA 0.00000000 1.54700725 1.48987379 2.891842 3.927415
image(bas.model, rotate=F)
coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0
## [1] TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE
confint(coef.model)
## 2.5% 97.5% beta
## Intercept 1.116139e+00 1.708031e+00 1.422090e+00
## avlength -8.564190e-02 1.467677e-02 -7.992245e-03
## avcondition -2.839253e+00 1.768409e+00 -9.325811e-02
## T_av -2.523565e-04 3.639709e-01 8.644408e-02
## O2_sat_av -2.287262e-02 1.596544e-02 -7.476618e-04
## Con_av -6.334651e-07 3.201875e-03 8.405554e-04
## COD_av 0.000000e+00 6.203673e-02 2.109188e-02
## NH4._av -4.940566e-01 1.660406e-01 -4.965441e-02
## Nt_av 0.000000e+00 2.277704e-01 5.507750e-02
## pool_riffle1 -4.507692e-04 9.323231e-01 1.389140e-01
## meander1 -5.556519e-01 2.496543e-01 -3.866036e-02
## netcen -2.956977e-05 2.000712e-05 -1.270941e-06
## updist -2.688456e-05 9.716094e-09 -5.952452e-06
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"
plot(confint(coef.model, parm = 2:11))
## NULL
confint <- confint(coef.model, parm = 2:11)
pip <- summary(bas.model)
PIP[c(1:12),12] <- pip[2:13,1]*sign(coef.model$postmean[2:13])
newdata1 <- newdata; newdata1[,"COD_av"] <- environment2$COD_av
BMA_IPI_COD_av <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3b = ggplot(environment2, aes(COD_av, BMA_IPI_COD_av$fit)) +
theme_bw() +
geom_line(color="red", size=1) +
geom_ribbon(aes(ymin = (BMA_IPI_COD_av$fit-BMA_IPI_COD_av$se.bma.fit), ymax = (BMA_IPI_COD_av$fit+BMA_IPI_COD_av$se.bma.fit)), alpha = .1) +
geom_point(data = environment2, aes(x=COD_av, y=avPI)) +
labs(x=expression("Chemical Oxygen Demand [mg/L]"), y=expression("I"[PI]*" - all parasites")) +
theme(axis.title.x = element_text(size=12),
axis.title.y = element_text(size=12)) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
labs(subtitle = "b)")
figure3b
bas.model <- bas.lm(avPI_ecto ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, prior="JZS")
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)
summary(bas.model)
## P(B != 0 | Y) model 1 model 2 model 3 model 4 model 5
## Intercept 1.0000000 1.0000000 1.000000000 1.0000000 1.000000 1.000000
## avlength 0.4669536 0.0000000 1.000000000 0.0000000 0.000000 1.000000
## avcondition 0.3758074 0.0000000 1.000000000 0.0000000 0.000000 0.000000
## T_av 0.2753561 0.0000000 1.000000000 0.0000000 0.000000 0.000000
## O2_sat_av 0.2919614 0.0000000 1.000000000 0.0000000 0.000000 0.000000
## Con_av 0.6617062 0.0000000 1.000000000 1.0000000 0.000000 0.000000
## COD_av 0.8437205 0.0000000 1.000000000 1.0000000 1.000000 1.000000
## NH4._av 0.4163705 0.0000000 1.000000000 0.0000000 0.000000 0.000000
## Nt_av 0.8302866 1.0000000 1.000000000 0.0000000 1.000000 1.000000
## pool_riffle1 0.4025108 0.0000000 1.000000000 0.0000000 0.000000 0.000000
## meander1 0.5031625 0.0000000 1.000000000 0.0000000 0.000000 0.000000
## netcen 0.3010717 0.0000000 1.000000000 0.0000000 0.000000 0.000000
## updist 0.3285744 0.0000000 1.000000000 0.0000000 0.000000 0.000000
## BF NA 0.1318791 0.008313608 0.4282678 0.319205 1.000000
## PostProbs NA 0.0510000 0.038600000 0.0301000 0.022500 0.021100
## R2 NA 0.3007000 0.682800000 0.4141000 0.403700 0.496600
## dim NA 2.0000000 13.000000000 3.0000000 3.000000 4.000000
## logmarg NA 4.0635523 1.299560343 5.2414154 4.947500 6.089422
image(bas.model, rotate=F)
coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0
## [1] TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE
confint(coef.model)
## 2.5% 97.5% beta
## Intercept 6.899850e-01 1.159346e+00 9.199577e-01
## avlength -1.072538e-01 6.599823e-03 -2.357999e-02
## avcondition -4.610896e+00 3.985984e-01 -7.055549e-01
## T_av -9.506040e-02 1.406152e-01 5.818763e-03
## O2_sat_av -1.638007e-02 2.228484e-02 1.122707e-03
## Con_av 0.000000e+00 3.193247e-03 1.203349e-03
## COD_av 0.000000e+00 6.428436e-02 3.405818e-02
## NH4._av -4.990487e-01 7.994328e-02 -8.034086e-02
## Nt_av -3.936480e-05 2.546442e-01 1.257830e-01
## pool_riffle1 -7.124275e-02 8.162781e-01 1.490513e-01
## meander1 -9.338322e-01 2.458840e-02 -2.358742e-01
## netcen -3.428514e-05 2.027454e-05 -2.447510e-06
## updist -1.553765e-05 4.040408e-06 -1.836166e-06
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"
plot(confint(coef.model, parm = 2:11))
## NULL
confint <- confint(coef.model, parm = 2:11)
pip <- summary(bas.model)
PIP[c(1:12),13] <- pip[2:13,1]*sign(coef.model$postmean[2:13])
newdata1 <- newdata; newdata1[,"COD_av"] <- environment2$COD_av
BMA_IPIecto_COD_av <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3c = ggplot(environment2, aes(COD_av, BMA_IPIecto_COD_av$fit)) +
theme_bw() +
geom_line(color="red", size=1) +
geom_ribbon(aes(ymin = (BMA_IPIecto_COD_av$fit-BMA_IPIecto_COD_av$se.bma.fit), ymax = (BMA_IPIecto_COD_av$fit+BMA_IPIecto_COD_av$se.bma.fit)), alpha = .1) +
geom_point(data = environment2, aes(x=COD_av, y=avPI_ecto)) +
labs(x=expression("Chemical Oxygen Demand [mg/L]"), y=expression("I"[PI]*" - ectoparasites")) +
theme(axis.title.x = element_text(size=12),
axis.title.y = element_text(size=12)) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
labs(subtitle = "c)")
figure3c
newdata1 <- newdata; newdata1[,"Nt_av"] <- environment2$Nt_av
BMA_IPIecto_Nt_av <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3f = ggplot(environment2, aes(Nt_av, BMA_IPIecto_Nt_av$fit)) +
theme_bw() +
geom_line(color="red", size=1) +
geom_ribbon(aes(ymin = (BMA_IPIecto_Nt_av$fit-BMA_IPIecto_Nt_av$se.bma.fit), ymax = (BMA_IPIecto_Nt_av$fit+BMA_IPIecto_Nt_av$se.bma.fit)), alpha = .1) +
geom_point(data = environment2, aes(x=Nt_av, y=avPI_ecto)) +
labs(x=expression("Nt [mg/L]"), y=expression("I"[PI]*" - ectoparasites")) +
theme(axis.title.x = element_text(size=12),
axis.title.y = element_text(size=12)) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
labs(subtitle = "f)")
figure3f
newdata1 <- newdata; newdata1[,"Con_av"] <- environment2$Con_av
BMA_IPIecto_Con_av <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3g = ggplot(environment2, aes(Con_av, BMA_IPIecto_Con_av$fit)) +
theme_bw() +
geom_line(color="red", size=1) +
geom_ribbon(aes(ymin = (BMA_IPIecto_Con_av$fit-BMA_IPIecto_Con_av$se.bma.fit), ymax = (BMA_IPIecto_Con_av$fit+BMA_IPIecto_Con_av$se.bma.fit)), alpha = .1) +
geom_point(data = environment2, aes(x=Con_av, y=avPI)) +
labs(x=expression(paste("Conductivity [", mu, "S/cm]")), y=expression("I"[PI]*" - all parasites")) +
theme(axis.title.x = element_text(size=12),
axis.title.y = element_text(size=12)) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
labs(subtitle = "g)")
figure3g
newdata1 <- newdata; newdata1[,"meander"] <- environment2$meander
BMA_IPIecto_meander <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3l = ggplot(environment2, aes(meander, BMA_IPIecto_meander$fit)) +
theme_bw() +
#geom_line(color="red", size=1) +
#geom_ribbon(aes(ymin = (BMA$fit-BMA$se.bma.fit), ymax = (BMA$fit+BMA$se.bma.fit)), alpha = .1) +
geom_point(data = environment2, aes(x=meander, y=avPI)) +
geom_boxplot(aes(lower = (BMA_IPIecto_meander$fit-BMA_IPIecto_meander$se.bma.fit), middle = BMA_IPIecto_meander$fit, upper = (BMA_IPIecto_meander$fit+BMA_IPIecto_meander$se.bma.fit))) +
labs(x=expression("Meander"), y=expression("I"[PI]*" - ectoparasites")) +
theme(axis.title.x = element_text(size=12),
axis.title.y = element_text(size=12)) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
labs(subtitle = "l)")
figure3l
bas.model <- bas.lm(avPI_endo ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, prior="JZS")
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)
summary(bas.model)
## P(B != 0 | Y) model 1 model 2 model 3 model 4 model 5
## Intercept 1.00000000 1.0000000 1.000000 1.00000000 1.0000000 1.0000000
## avlength 0.06442140 0.0000000 0.000000 0.00000000 0.0000000 0.0000000
## avcondition 0.12246410 0.0000000 0.000000 1.00000000 0.0000000 0.0000000
## T_av 0.30397004 0.0000000 1.000000 0.00000000 0.0000000 1.0000000
## O2_sat_av 0.05178919 0.0000000 0.000000 0.00000000 0.0000000 0.0000000
## Con_av 0.05158494 0.0000000 0.000000 0.00000000 0.0000000 0.0000000
## COD_av 0.04593088 0.0000000 0.000000 0.00000000 0.0000000 0.0000000
## NH4._av 0.08528364 0.0000000 0.000000 0.00000000 0.0000000 1.0000000
## Nt_av 0.06477759 0.0000000 0.000000 0.00000000 0.0000000 0.0000000
## pool_riffle1 0.04504579 0.0000000 0.000000 0.00000000 0.0000000 0.0000000
## meander1 0.05899866 0.0000000 0.000000 0.00000000 0.0000000 0.0000000
## netcen 0.05414430 0.0000000 0.000000 0.00000000 0.0000000 0.0000000
## updist 0.10936917 0.0000000 0.000000 0.00000000 1.0000000 0.0000000
## BF NA 0.3484238 1.000000 0.33059983 0.2877395 0.6134470
## PostProbs NA 0.4921000 0.117700 0.03890000 0.0339000 0.0131000
## R2 NA 0.0000000 0.166100 0.11040000 0.1032000 0.2244000
## dim NA 1.0000000 2.000000 2.00000000 2.0000000 3.0000000
## logmarg NA 0.0000000 1.054336 -0.05251095 -0.1913639 0.5656742
image(bas.model, rotate=F)
coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0
## [1] TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE
confint(coef.model)
## 2.5% 97.5% beta
## Intercept 2.974020e-01 6.917548e-01 4.947994e-01
## avlength -3.456166e-04 1.669276e-02 1.398963e-03
## avcondition -3.455674e-03 2.596667e+00 2.446911e-01
## T_av 0.000000e+00 2.137081e-01 4.520999e-02
## O2_sat_av -3.702737e-04 3.813693e-05 -2.446396e-04
## Con_av -6.714701e-05 4.289287e-06 -1.385311e-05
## COD_av 0.000000e+00 0.000000e+00 -1.072193e-04
## NH4._av -1.119913e-01 1.101911e-04 -9.322620e-03
## Nt_av -2.836973e-02 1.844054e-03 -2.149910e-03
## pool_riffle1 0.000000e+00 0.000000e+00 2.267863e-03
## meander1 -1.340521e-02 7.816289e-02 1.005422e-02
## netcen -1.286967e-06 8.497622e-08 -3.607730e-07
## updist -9.420533e-06 0.000000e+00 -8.091451e-07
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"
plot(confint(coef.model, parm = 2:11))
## Warning in arrows(x[not.deg], ci[not.deg, 1], x[not.deg], ci[not.deg, 2], :
## zero-length arrow is of indeterminate angle and so skipped
## NULL
confint <- confint(coef.model, parm = 2:11)
pip <- summary(bas.model)
PIP[c(1:12),14] <- pip[2:13,1]*sign(coef.model$postmean[2:13])
library(gridExtra)
pdf("Figure3.pdf", width = 10.8, height = 14.4)
grid.arrange(figure3a, figure3b, figure3c, figure3d, figure3e, figure3f, figure3g, figure3h, figure3i, figure3j, figure3k, figure3l)
dev.off()
## png
## 2
Model-based analysis of multivariate abundance data using Bayesian Markov chain Monte Carlo methods for parameter estimation
data$Site <- as.factor(data$site)
levels(data$site) <- levels(as.factor(environment2$Site))
data_m <- merge(data, environment2, by = "Site")
data_all <- na.omit(data_m)
names(data_all)
## [1] "Site" "site" "fish"
## [4] "parspeciesrichness" "div_shannon" "div_simpson"
## [7] "temp" "pH" "conductivity"
## [10] "nitrogen" "phosphorus" "oxygen"
## [13] "netcen.x" "updist.x" "updist2"
## [16] "updist3" "fishspeciesrichness" "weight"
## [19] "weigh..g." "length" "SMI"
## [22] "Sex" "Gyr" "Tri"
## [25] "Glu" "ecto_screener" "Con"
## [28] "CysL" "Pro" "Aca"
## [31] "Cam" "Ang" "CysI"
## [34] "endo_screener" "PI" "PI_ecto"
## [37] "PI_endo" "T_av" "O2_sat_av"
## [40] "Con_av" "COD_av" "NH4._av"
## [43] "Nt_av" "SM_av" "pool_riffle"
## [46] "meander" "updist.y" "netcen.y"
avcondition <- aggregate(data$SMI, by = list(data[,1]), function(x){mean(x, na.rm =T)})[,2]
avlength <- aggregate(data$length, by = list(data[,1]), function(x){mean(x, na.rm =T)})[,2]
y <- round(cbind(avab$Gyr, avab$Tri, avab$Glu, avab$Con, avab$Ang))
X <- cbind(avcondition,
avlength,
environment2$T_av,
environment2$O2_sat_av,
environment2$Con_av,
environment2$COD_av,
environment2$NH4._av,
environment2$Nt_av,
environment2$netcen,
environment2$updist,
as.numeric(environment2$pool_riffle),
as.numeric(environment2$meander))
colnames(X) <- c("avcondition", "avlength", "T", "O2", "Con", "COD", "NH4", "Nt", "netcen", "updist", "pool_riffle", "meander")
example_mcmc_control <- list(n.burnin = 1000, n.iteration = 10000, n.thin = 1)
testpath <- file.path(tempdir(), "jagsboralmodel.txt")
paramod <- boral(y, X = X,
family = "negative.binomial",
mcmc.control = example_mcmc_control,
model.name = testpath,
lv.control = list(num.lv = 2, type = "independent"),
save.model = TRUE)
## module glm loaded
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 185
## Unobserved stochastic nodes: 338
## Total graph size: 2173
##
## Initializing model
plot(paramod)
## NULL
coefsplot(covname = "avcondition", object = paramod) #Condition
coefsplot(covname = "avlength", object = paramod) #Length
coefsplot(covname = "T", object = paramod) #Temperature
coefsplot(covname = "O2", object = paramod) #Oxygen
coefsplot(covname = "Con", object = paramod) #Conductivity
coefsplot(covname = "COD", object = paramod) #COD
coefsplot(covname = "NH4", object = paramod) #NH4
coefsplot(covname = "Nt", object = paramod) #Nt
coefsplot(covname = "netcen", object = paramod) #netcen
coefsplot(covname = "updist", object = paramod) #updist
coefsplot(covname = "pool_riffle", object = paramod) #poolriffle
coefsplot(covname = "meander", object = paramod) #meander
envcors <- get.enviro.cor(paramod)
rescors <- get.residual.cor(paramod)
corrplot(envcors$sig.cor, type = "lower", diag = FALSE, title = "Correlations due to covariates", mar = c(3,0.5,2,1), tl.srt = 45)
corrplot(rescors$sig.cor, type = "lower", diag = FALSE, title = "Residual correlations", mar = c(3,0.5,2,1), tl.srt = 45)
y <- round(cbind(medin$Gyr, medin$Tri, medin$Glu, medin$Con, medin$Ang))
paramod <- boral(y, X = X,
family = "negative.binomial",
mcmc.control = example_mcmc_control,
model.name = testpath,
lv.control = list(num.lv = 2, type = "independent"),
save.model = TRUE)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 185
## Unobserved stochastic nodes: 338
## Total graph size: 2173
##
## Initializing model
plot(paramod)
## NULL
coefsplot(covname = "avcondition", object = paramod) #Condition
coefsplot(covname = "avlength", object = paramod) #Length
coefsplot(covname = "T", object = paramod) #Temperature
coefsplot(covname = "O2", object = paramod) #Oxygen
coefsplot(covname = "Con", object = paramod) #Conductivity
coefsplot(covname = "COD", object = paramod) #COD
coefsplot(covname = "NH4", object = paramod) #NH4
coefsplot(covname = "Nt", object = paramod) #Nt
coefsplot(covname = "netcen", object = paramod) #netcen
coefsplot(covname = "updist", object = paramod) #updist
coefsplot(covname = "pool_riffle", object = paramod) #poolriffle
coefsplot(covname = "meander", object = paramod) #meander
envcors <- get.enviro.cor(paramod)
rescors <- get.residual.cor(paramod)
corrplot(envcors$sig.cor, type = "lower", diag = FALSE, title = "Correlations due to covariates", mar = c(3,0.5,2,1), tl.srt = 45)
corrplot(rescors$sig.cor, type = "lower", diag = FALSE, title = "Residual correlations", mar = c(3,0.5,2,1), tl.srt = 45)
# Component communities: Bray-Curtis dissimilarities based on Hellinger transformed average abundance data
spe.hel_bray <- vegdist(decostand(avab[,-1], na.rm=T, method="hellinger"), method="bray", na.rm=T)
# Check whether Euclidean and Bray-Curtis distances are comparable
spe.hel_euc <- vegdist(decostand(avab[,-1], na.rm=T, method="hellinger"), method="euc", na.rm=T)
plot(spe.hel_bray, spe.hel_euc)
mantel(spe.hel_bray, spe.hel_euc)
##
## Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel(xdis = spe.hel_bray, ydis = spe.hel_euc)
##
## Mantel statistic r: 0.9648
## Significance: 0.001
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.105 0.138 0.165 0.196
## Permutation: free
## Number of permutations: 999
model_adonis = adonis(spe.hel_bray ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2)
## 'adonis' will be deprecated: use 'adonis2' instead
model_adonis$aov.tab
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## avlength 1 0.1863 0.18629 1.8398 0.04303 0.117
## avcondition 1 0.1305 0.13055 1.2893 0.03016 0.267
## T_av 1 0.3204 0.32039 3.1643 0.07401 0.016 *
## O2_sat_av 1 0.1368 0.13678 1.3509 0.03160 0.241
## Con_av 1 0.1326 0.13262 1.3098 0.03064 0.266
## COD_av 1 0.0657 0.06568 0.6486 0.01517 0.668
## NH4._av 1 0.2040 0.20399 2.0146 0.04712 0.096 .
## Nt_av 1 0.1451 0.14509 1.4329 0.03352 0.227
## pool_riffle 1 0.0853 0.08529 0.8424 0.01970 0.517
## meander 1 0.2029 0.20286 2.0035 0.04686 0.098 .
## netcen 1 0.2173 0.21728 2.1459 0.05019 0.088 .
## updist 1 0.0719 0.07193 0.7104 0.01662 0.621
## Residuals 24 2.4301 0.10125 0.56137
## Total 36 4.3288 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# environmental variables
env_select <- environment2[,c("T_av", "O2_sat_av", "Con_av", "COD_av", "NH4._av", "Nt_av", "pool_riffle", "meander", "netcen", "updist")]
env_select$pool_riffle <- as.numeric(env_select$pool_riffle)
env_select$meander <- as.numeric(env_select$meander)
pca <- prcomp(env_select, scale.=T)
summary(pca)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 1.7124 1.5545 1.1221 1.0140 0.88807 0.79463 0.56647
## Proportion of Variance 0.2933 0.2416 0.1259 0.1028 0.07887 0.06314 0.03209
## Cumulative Proportion 0.2933 0.5349 0.6608 0.7636 0.84248 0.90563 0.93771
## PC8 PC9 PC10
## Standard deviation 0.50483 0.46939 0.38429
## Proportion of Variance 0.02549 0.02203 0.01477
## Cumulative Proportion 0.96320 0.98523 1.00000
plot(pca)
biplot(pca)
# Assess the effect of environmental variables on parasite component community dissimilarities using distance based RDA
spe.rda <- dbrda(spe.hel_bray ~ T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander, data = env_select)
summary(spe.rda)
##
## Call:
## dbrda(formula = spe.hel_bray ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander, data = env_select)
##
## Partitioning of squared Bray distance:
## Inertia Proportion
## Total 4.329 1.0000
## Constrained 1.166 0.2695
## Unconstrained 3.162 0.7305
##
## Eigenvalues, and their contribution to the squared Bray distance
##
## Importance of components:
## dbRDA1 dbRDA2 dbRDA3 dbRDA4 dbRDA5 dbRDA6
## Eigenvalue 0.5740 0.3372 0.17566 0.09415 0.042525 0.0020670
## Proportion Explained 0.1326 0.0779 0.04058 0.02175 0.009824 0.0004775
## Cumulative Proportion NA NA NA NA NA NA
## idbRDA1 idbRDA2 MDS1 MDS2 MDS3 MDS4 MDS5
## Eigenvalue -0.014797 -0.04441 1.3249 0.8027 0.4680 0.31323 0.29121
## Proportion Explained 0.003418 0.01026 0.3061 0.1854 0.1081 0.07236 0.06727
## Cumulative Proportion NA NA NA NA NA NA NA
## MDS6 MDS7 MDS8 MDS9 MDS10 MDS11 MDS12
## Eigenvalue 0.13644 0.11438 0.08987 0.07065 0.06425 0.02766 0.01576
## Proportion Explained 0.03152 0.02642 0.02076 0.01632 0.01484 0.00639 0.00364
## Cumulative Proportion NA NA NA NA NA NA NA
## MDS13 MDS14 iMDS1 iMDS2 iMDS3
## Eigenvalue 0.011872 2.725e-04 -0.0020314 -0.007356 -0.01078
## Proportion Explained 0.002743 6.295e-05 0.0004693 0.001699 0.00249
## Cumulative Proportion NA NA NA NA NA
## iMDS4 iMDS5 iMDS6 iMDS7 iMDS8 iMDS9
## Eigenvalue -0.014744 -0.01935 -0.022495 -0.029849 -0.03342 -0.041489
## Proportion Explained 0.003406 0.00447 0.005197 0.006895 0.00772 0.009584
## Cumulative Proportion NA NA NA NA NA NA
## iMDS10 iMDS11 iMDS12 iMDS13 iMDS14
## Eigenvalue -0.05984 -0.06439 -0.06995 -0.09116 -0.10203
## Proportion Explained 0.01382 0.01487 0.01616 0.02106 0.02357
## Cumulative Proportion NA NA NA NA NA
##
## Accumulated constrained eigenvalues
## Importance of components:
## dbRDA1 dbRDA2 dbRDA3 dbRDA4 dbRDA5 dbRDA6 idbRDA1
## Eigenvalue 0.5740 0.3372 0.1757 0.09415 0.04252 0.002067 -0.01480
## Proportion Explained 0.4921 0.2891 0.1506 0.08072 0.03646 0.001772 0.01269
## Cumulative Proportion NA NA NA NA NA NA NA
## idbRDA2
## Eigenvalue -0.04441
## Proportion Explained 0.03807
## Cumulative Proportion NA
##
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores: 3.533199
##
##
## Site scores (weighted sums of species scores)
##
## dbRDA1 dbRDA2 dbRDA3 dbRDA4 dbRDA5 dbRDA6
## 1 -0.07530 -0.62580 -0.63602 -2.643196 -0.34427 -23.8946
## 2 0.87242 0.77832 -0.32068 2.074749 -0.25800 41.3754
## 3 -0.81094 1.29200 -0.57735 -0.708420 -2.35800 -5.0035
## 4 -0.97450 0.91051 1.79849 -2.192492 -1.99248 -4.7134
## 5 -0.25889 0.74559 -1.30995 -1.418227 0.33185 -21.9720
## 6 0.41640 1.01888 -0.34551 -2.318324 1.61924 -24.8949
## 7 -0.98286 -0.67495 0.62017 1.174218 1.55155 -0.8974
## 8 0.93807 0.62635 -3.10878 -0.765073 2.41981 -3.0764
## 9 -0.62816 -0.77245 -1.30540 1.806170 1.30543 -41.2679
## 10 -1.03727 -0.50089 0.62572 0.369727 1.00378 4.7819
## 11 1.10743 -0.99448 -0.66515 -0.071830 0.83402 -14.8528
## 12 0.61624 0.76828 0.87110 0.357836 0.12125 38.2552
## 13 -0.96088 -0.72453 0.47977 1.569869 1.43399 -6.8692
## 14 1.41128 2.03453 -1.67256 0.776255 0.27648 29.9435
## 15 -1.00448 -0.59302 0.65010 0.833950 1.44643 3.4902
## 16 1.26382 0.81623 3.19230 1.380128 1.12523 -10.4959
## 17 -0.92726 -0.73637 0.33774 1.892276 1.07787 -8.2400
## 18 -1.04523 -1.01280 -0.59052 2.897871 0.35870 -14.2865
## 19 0.89987 -1.54255 -0.54435 1.738188 -2.61939 32.8867
## 20 -0.14608 -0.34079 0.78630 1.316880 1.70249 -9.1246
## 21 -1.30124 -1.46318 -1.27370 3.739798 0.91928 -16.7232
## 22 -0.07485 1.73636 0.30728 -1.760677 -0.46041 6.1504
## 23 -0.18008 -0.64105 0.55428 -1.310793 -0.82437 -3.6579
## 24 0.53054 0.08088 1.31623 -3.529128 -1.96164 41.6937
## 25 -0.50901 -0.82222 -1.26480 2.167458 -6.71178 -22.4535
## 26 0.54874 0.92050 1.12409 0.009258 0.07714 34.0773
## 27 1.25804 -0.69768 2.37255 -1.697654 -3.24936 -16.2349
## 28 -0.86554 0.32492 -0.68347 0.704689 -0.16329 13.0277
## 29 1.82276 0.98120 0.33891 -1.852448 4.33581 -8.6587
## 30 0.72974 -0.59907 -1.00943 -0.201449 1.77124 -4.5186
## 31 -1.02904 0.19515 1.05102 -1.468951 -1.36610 23.1359
## 32 -0.13877 -0.06050 -0.19569 -1.899083 -1.10403 -31.4856
## 33 1.01209 -0.59332 2.58400 -0.318337 -1.43566 39.0635
## 34 -0.39319 -1.79743 0.20497 -0.362972 -0.59693 -19.8004
## 35 -0.89506 1.46497 -0.05622 -1.876722 -1.86120 3.5844
## 36 0.87386 1.17335 -3.50567 1.196922 1.51036 -9.4970
## 37 -0.06267 -0.67495 -0.14975 0.389535 2.08495 11.1529
##
##
## Site constraints (linear combinations of constraining variables)
##
## dbRDA1 dbRDA2 dbRDA3 dbRDA4 dbRDA5 dbRDA6
## 1 -0.16605 -0.69736 -0.45052 0.39949 -1.140870 0.26766
## 2 -0.07683 -0.33692 0.61638 0.50606 -0.381728 1.66676
## 3 -0.10071 0.76299 0.51759 0.02910 0.003366 -0.26259
## 4 -1.55548 0.90013 0.66933 -0.51241 -0.415318 -0.13093
## 5 -0.41914 0.27994 -0.66104 -1.10909 0.384267 0.50060
## 6 0.63090 0.30135 -0.07724 0.37780 0.857007 -0.33310
## 7 -0.15252 -0.25647 -0.20458 -0.38652 0.792931 0.48033
## 8 -0.14874 0.12938 -0.39866 -0.78891 0.083814 0.27677
## 9 -0.14359 0.11628 0.95439 -0.04629 0.609992 -0.64671
## 10 -0.74154 -0.65088 0.87072 -1.17151 -0.875610 -1.73263
## 11 0.77571 -0.48314 0.02321 -0.76003 0.259283 -0.01494
## 12 -0.44387 -0.28442 0.54938 -0.01813 1.062736 -0.71676
## 13 -0.25608 -0.61494 -0.47847 0.74269 0.180828 -0.14346
## 14 -0.23636 0.67472 -0.87185 0.72759 0.046993 -0.61306
## 15 -0.75508 0.19781 0.41035 0.89449 -0.744271 -0.25040
## 16 0.97946 1.50039 0.74911 0.12808 -0.657115 -0.30363
## 17 -1.19679 -0.18091 -0.90174 0.16093 -0.229189 0.70454
## 18 0.23197 0.39183 -0.42932 -0.82808 -0.242760 0.64991
## 19 0.80868 -0.55920 -0.22435 0.32426 -0.445610 -0.51746
## 20 -0.19206 0.12003 0.24728 0.99541 0.747110 0.40576
## 21 -0.38477 0.09371 -0.03632 1.07969 1.383306 -0.25819
## 22 0.02352 0.34014 0.07462 0.03662 -0.352806 0.35626
## 23 0.18716 -0.91860 -0.02126 -0.69548 0.117329 0.60753
## 24 0.03851 0.57991 -0.70662 -1.10691 0.533700 -0.10095
## 25 0.61401 -0.01017 -0.51875 0.25795 -1.076732 -0.24901
## 26 0.42275 0.07490 0.99168 0.37950 -0.249480 0.19511
## 27 0.54478 -0.12356 0.31406 0.41569 0.237094 0.09772
## 28 0.02044 0.43106 -0.71781 0.33274 -0.428303 -0.26278
## 29 1.23235 0.15293 0.10463 -0.07285 0.107338 0.07885
## 30 0.02214 -1.33353 -0.37271 0.38762 -0.062960 -0.48751
## 31 -0.71637 0.11204 -0.85950 0.24734 -0.066602 -0.17455
## 32 0.01756 -0.71899 -0.27770 0.01727 0.106473 -0.73172
## 33 0.78296 -0.75806 0.27650 -0.27767 -0.031297 -0.01895
## 34 -0.33590 -0.63554 1.07989 -0.07860 -0.526645 0.97933
## 35 -0.15761 0.49842 0.66426 -0.05091 0.331527 0.71522
## 36 0.53360 0.88375 -0.84613 0.03571 -0.630019 -0.08456
## 37 0.31300 0.02099 -0.05877 -0.57261 0.712218 0.05155
##
##
## Biplot scores for constraining variables
##
## dbRDA1 dbRDA2 dbRDA3 dbRDA4 dbRDA5 dbRDA6
## T_av 0.63982 0.25527 0.1662 -0.34674 -0.2289 -0.48193
## O2_sat_av 0.02371 -0.02938 -0.8854 -0.18557 0.2035 -0.13356
## Con_av -0.25315 -0.27391 -0.1407 -0.76998 -0.4264 0.07657
## COD_av -0.13235 0.15034 0.4807 0.08407 0.1430 -0.57341
## NH4._av -0.40051 0.14159 0.7569 -0.21438 -0.2051 -0.37005
## Nt_av -0.31714 -0.23687 0.2952 -0.33830 -0.5347 -0.40917
## pool_riffle -0.34609 -0.15029 -0.4384 0.27287 -0.5547 -0.42271
## meander -0.33563 0.52703 -0.3800 -0.21165 -0.5422 -0.13159
anova(spe.rda)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: dbrda(formula = spe.hel_bray ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander, data = env_select)
## Df SumOfSqs F Pr(>F)
## Model 8 1.1664 1.2909 0.18
## Residual 28 3.1624
RsquareAdj(spe.rda)$adj.r.squared
## [1] 0.06072755
mod0 <- dbrda(spe.hel_bray ~ 1, env_select[,-c(9:10)]) # Model with intercept only #edit_PH
mod1 <- dbrda(spe.hel_bray ~ ., env_select[,-c(9:10)]) # Model with all explanatory variables #edit_PH
step.res <- ordiR2step(mod0, mod1, direction = "both",perm.max = 200)
## Step: R2.adj= 0
## Call: spe.hel_bray ~ 1
##
## R2.adjusted
## <All variables> 0.060727547
## + T_av 0.036398792
## + NH4._av 0.020208612
## + meander 0.018502880
## + O2_sat_av 0.004277611
## + Con_av 0.001872668
## + pool_riffle 0.001742860
## <none> 0.000000000
## + Nt_av -0.002060170
## + COD_av -0.017968936
##
## Df AIC F Pr(>F)
## + T_av 1 54.788 2.3599 0.052 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
step.res$anova # Summary table
## NULL
plot(spe.rda, scaling = 1) # it is for technical reasons not possible to plot both site and species scores
summary(spe.rda)
##
## Call:
## dbrda(formula = spe.hel_bray ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander, data = env_select)
##
## Partitioning of squared Bray distance:
## Inertia Proportion
## Total 4.329 1.0000
## Constrained 1.166 0.2695
## Unconstrained 3.162 0.7305
##
## Eigenvalues, and their contribution to the squared Bray distance
##
## Importance of components:
## dbRDA1 dbRDA2 dbRDA3 dbRDA4 dbRDA5 dbRDA6
## Eigenvalue 0.5740 0.3372 0.17566 0.09415 0.042525 0.0020670
## Proportion Explained 0.1326 0.0779 0.04058 0.02175 0.009824 0.0004775
## Cumulative Proportion NA NA NA NA NA NA
## idbRDA1 idbRDA2 MDS1 MDS2 MDS3 MDS4 MDS5
## Eigenvalue -0.014797 -0.04441 1.3249 0.8027 0.4680 0.31323 0.29121
## Proportion Explained 0.003418 0.01026 0.3061 0.1854 0.1081 0.07236 0.06727
## Cumulative Proportion NA NA NA NA NA NA NA
## MDS6 MDS7 MDS8 MDS9 MDS10 MDS11 MDS12
## Eigenvalue 0.13644 0.11438 0.08987 0.07065 0.06425 0.02766 0.01576
## Proportion Explained 0.03152 0.02642 0.02076 0.01632 0.01484 0.00639 0.00364
## Cumulative Proportion NA NA NA NA NA NA NA
## MDS13 MDS14 iMDS1 iMDS2 iMDS3
## Eigenvalue 0.011872 2.725e-04 -0.0020314 -0.007356 -0.01078
## Proportion Explained 0.002743 6.295e-05 0.0004693 0.001699 0.00249
## Cumulative Proportion NA NA NA NA NA
## iMDS4 iMDS5 iMDS6 iMDS7 iMDS8 iMDS9
## Eigenvalue -0.014744 -0.01935 -0.022495 -0.029849 -0.03342 -0.041489
## Proportion Explained 0.003406 0.00447 0.005197 0.006895 0.00772 0.009584
## Cumulative Proportion NA NA NA NA NA NA
## iMDS10 iMDS11 iMDS12 iMDS13 iMDS14
## Eigenvalue -0.05984 -0.06439 -0.06995 -0.09116 -0.10203
## Proportion Explained 0.01382 0.01487 0.01616 0.02106 0.02357
## Cumulative Proportion NA NA NA NA NA
##
## Accumulated constrained eigenvalues
## Importance of components:
## dbRDA1 dbRDA2 dbRDA3 dbRDA4 dbRDA5 dbRDA6 idbRDA1
## Eigenvalue 0.5740 0.3372 0.1757 0.09415 0.04252 0.002067 -0.01480
## Proportion Explained 0.4921 0.2891 0.1506 0.08072 0.03646 0.001772 0.01269
## Cumulative Proportion NA NA NA NA NA NA NA
## idbRDA2
## Eigenvalue -0.04441
## Proportion Explained 0.03807
## Cumulative Proportion NA
##
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores: 3.533199
##
##
## Site scores (weighted sums of species scores)
##
## dbRDA1 dbRDA2 dbRDA3 dbRDA4 dbRDA5 dbRDA6
## 1 -0.07530 -0.62580 -0.63602 -2.643196 -0.34427 -23.8946
## 2 0.87242 0.77832 -0.32068 2.074749 -0.25800 41.3754
## 3 -0.81094 1.29200 -0.57735 -0.708420 -2.35800 -5.0035
## 4 -0.97450 0.91051 1.79849 -2.192492 -1.99248 -4.7134
## 5 -0.25889 0.74559 -1.30995 -1.418227 0.33185 -21.9720
## 6 0.41640 1.01888 -0.34551 -2.318324 1.61924 -24.8949
## 7 -0.98286 -0.67495 0.62017 1.174218 1.55155 -0.8974
## 8 0.93807 0.62635 -3.10878 -0.765073 2.41981 -3.0764
## 9 -0.62816 -0.77245 -1.30540 1.806170 1.30543 -41.2679
## 10 -1.03727 -0.50089 0.62572 0.369727 1.00378 4.7819
## 11 1.10743 -0.99448 -0.66515 -0.071830 0.83402 -14.8528
## 12 0.61624 0.76828 0.87110 0.357836 0.12125 38.2552
## 13 -0.96088 -0.72453 0.47977 1.569869 1.43399 -6.8692
## 14 1.41128 2.03453 -1.67256 0.776255 0.27648 29.9435
## 15 -1.00448 -0.59302 0.65010 0.833950 1.44643 3.4902
## 16 1.26382 0.81623 3.19230 1.380128 1.12523 -10.4959
## 17 -0.92726 -0.73637 0.33774 1.892276 1.07787 -8.2400
## 18 -1.04523 -1.01280 -0.59052 2.897871 0.35870 -14.2865
## 19 0.89987 -1.54255 -0.54435 1.738188 -2.61939 32.8867
## 20 -0.14608 -0.34079 0.78630 1.316880 1.70249 -9.1246
## 21 -1.30124 -1.46318 -1.27370 3.739798 0.91928 -16.7232
## 22 -0.07485 1.73636 0.30728 -1.760677 -0.46041 6.1504
## 23 -0.18008 -0.64105 0.55428 -1.310793 -0.82437 -3.6579
## 24 0.53054 0.08088 1.31623 -3.529128 -1.96164 41.6937
## 25 -0.50901 -0.82222 -1.26480 2.167458 -6.71178 -22.4535
## 26 0.54874 0.92050 1.12409 0.009258 0.07714 34.0773
## 27 1.25804 -0.69768 2.37255 -1.697654 -3.24936 -16.2349
## 28 -0.86554 0.32492 -0.68347 0.704689 -0.16329 13.0277
## 29 1.82276 0.98120 0.33891 -1.852448 4.33581 -8.6587
## 30 0.72974 -0.59907 -1.00943 -0.201449 1.77124 -4.5186
## 31 -1.02904 0.19515 1.05102 -1.468951 -1.36610 23.1359
## 32 -0.13877 -0.06050 -0.19569 -1.899083 -1.10403 -31.4856
## 33 1.01209 -0.59332 2.58400 -0.318337 -1.43566 39.0635
## 34 -0.39319 -1.79743 0.20497 -0.362972 -0.59693 -19.8004
## 35 -0.89506 1.46497 -0.05622 -1.876722 -1.86120 3.5844
## 36 0.87386 1.17335 -3.50567 1.196922 1.51036 -9.4970
## 37 -0.06267 -0.67495 -0.14975 0.389535 2.08495 11.1529
##
##
## Site constraints (linear combinations of constraining variables)
##
## dbRDA1 dbRDA2 dbRDA3 dbRDA4 dbRDA5 dbRDA6
## 1 -0.16605 -0.69736 -0.45052 0.39949 -1.140870 0.26766
## 2 -0.07683 -0.33692 0.61638 0.50606 -0.381728 1.66676
## 3 -0.10071 0.76299 0.51759 0.02910 0.003366 -0.26259
## 4 -1.55548 0.90013 0.66933 -0.51241 -0.415318 -0.13093
## 5 -0.41914 0.27994 -0.66104 -1.10909 0.384267 0.50060
## 6 0.63090 0.30135 -0.07724 0.37780 0.857007 -0.33310
## 7 -0.15252 -0.25647 -0.20458 -0.38652 0.792931 0.48033
## 8 -0.14874 0.12938 -0.39866 -0.78891 0.083814 0.27677
## 9 -0.14359 0.11628 0.95439 -0.04629 0.609992 -0.64671
## 10 -0.74154 -0.65088 0.87072 -1.17151 -0.875610 -1.73263
## 11 0.77571 -0.48314 0.02321 -0.76003 0.259283 -0.01494
## 12 -0.44387 -0.28442 0.54938 -0.01813 1.062736 -0.71676
## 13 -0.25608 -0.61494 -0.47847 0.74269 0.180828 -0.14346
## 14 -0.23636 0.67472 -0.87185 0.72759 0.046993 -0.61306
## 15 -0.75508 0.19781 0.41035 0.89449 -0.744271 -0.25040
## 16 0.97946 1.50039 0.74911 0.12808 -0.657115 -0.30363
## 17 -1.19679 -0.18091 -0.90174 0.16093 -0.229189 0.70454
## 18 0.23197 0.39183 -0.42932 -0.82808 -0.242760 0.64991
## 19 0.80868 -0.55920 -0.22435 0.32426 -0.445610 -0.51746
## 20 -0.19206 0.12003 0.24728 0.99541 0.747110 0.40576
## 21 -0.38477 0.09371 -0.03632 1.07969 1.383306 -0.25819
## 22 0.02352 0.34014 0.07462 0.03662 -0.352806 0.35626
## 23 0.18716 -0.91860 -0.02126 -0.69548 0.117329 0.60753
## 24 0.03851 0.57991 -0.70662 -1.10691 0.533700 -0.10095
## 25 0.61401 -0.01017 -0.51875 0.25795 -1.076732 -0.24901
## 26 0.42275 0.07490 0.99168 0.37950 -0.249480 0.19511
## 27 0.54478 -0.12356 0.31406 0.41569 0.237094 0.09772
## 28 0.02044 0.43106 -0.71781 0.33274 -0.428303 -0.26278
## 29 1.23235 0.15293 0.10463 -0.07285 0.107338 0.07885
## 30 0.02214 -1.33353 -0.37271 0.38762 -0.062960 -0.48751
## 31 -0.71637 0.11204 -0.85950 0.24734 -0.066602 -0.17455
## 32 0.01756 -0.71899 -0.27770 0.01727 0.106473 -0.73172
## 33 0.78296 -0.75806 0.27650 -0.27767 -0.031297 -0.01895
## 34 -0.33590 -0.63554 1.07989 -0.07860 -0.526645 0.97933
## 35 -0.15761 0.49842 0.66426 -0.05091 0.331527 0.71522
## 36 0.53360 0.88375 -0.84613 0.03571 -0.630019 -0.08456
## 37 0.31300 0.02099 -0.05877 -0.57261 0.712218 0.05155
##
##
## Biplot scores for constraining variables
##
## dbRDA1 dbRDA2 dbRDA3 dbRDA4 dbRDA5 dbRDA6
## T_av 0.63982 0.25527 0.1662 -0.34674 -0.2289 -0.48193
## O2_sat_av 0.02371 -0.02938 -0.8854 -0.18557 0.2035 -0.13356
## Con_av -0.25315 -0.27391 -0.1407 -0.76998 -0.4264 0.07657
## COD_av -0.13235 0.15034 0.4807 0.08407 0.1430 -0.57341
## NH4._av -0.40051 0.14159 0.7569 -0.21438 -0.2051 -0.37005
## Nt_av -0.31714 -0.23687 0.2952 -0.33830 -0.5347 -0.40917
## pool_riffle -0.34609 -0.15029 -0.4384 0.27287 -0.5547 -0.42271
## meander -0.33563 0.52703 -0.3800 -0.21165 -0.5422 -0.13159
anova(spe.rda)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: dbrda(formula = spe.hel_bray ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander, data = env_select)
## Df SumOfSqs F Pr(>F)
## Model 8 1.1664 1.2909 0.181
## Residual 28 3.1624
anova(spe.rda, by="term")
## Permutation test for dbrda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## Model: dbrda(formula = spe.hel_bray ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander, data = env_select)
## Df SumOfSqs F Pr(>F)
## T_av 1 0.2734 2.4210 0.061 .
## O2_sat_av 1 0.1377 1.2195 0.319
## Con_av 1 0.1613 1.4283 0.246
## COD_av 1 0.0676 0.5990 0.665
## NH4._av 1 0.1709 1.5131 0.229
## Nt_av 1 0.0501 0.4439 0.821
## pool_riffle 1 0.0657 0.5818 0.692
## meander 1 0.2396 2.1210 0.083 .
## Residual 28 3.1624
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova.cca(spe.rda, step=1000);
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: dbrda(formula = spe.hel_bray ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander, data = env_select)
## Df SumOfSqs F Pr(>F)
## Model 8 1.1664 1.2909 0.153
## Residual 28 3.1624
anova.cca(spe.rda, step=1000, by="term");
## Permutation test for dbrda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## Model: dbrda(formula = spe.hel_bray ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander, data = env_select)
## Df SumOfSqs F Pr(>F)
## T_av 1 0.2734 2.4210 0.063 .
## O2_sat_av 1 0.1377 1.2195 0.287
## Con_av 1 0.1613 1.4283 0.207
## COD_av 1 0.0676 0.5990 0.677
## NH4._av 1 0.1709 1.5131 0.212
## Nt_av 1 0.0501 0.4439 0.792
## pool_riffle 1 0.0657 0.5818 0.687
## meander 1 0.2396 2.1210 0.087 .
## Residual 28 3.1624
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RsquareAdj(spe.rda)$adj.r.squared;
## [1] 0.06072755
RsquareAdj(spe.rda)$r.squared
## [1] 0.2694548
# Generate a plot of the RDA
rownames(spe.rda$CCA$biplot) = c("temperature","O2","conductivity","COD","NH4","Nt","pool riffle pattern","meander")
ordiplot(spe.rda, choices=c(1,2), type="n", xlab="dbRDA 1 (49.21 %)", ylab="dbRDA 2 (28.91 %)")
points(spe.rda, "sites", pch=21, col="black", bg="white")
text(spe.rda, "cn", col="grey", cex=0.9, scaling=1)
# Save the plot to an object
figure4c = recordPlot()
# Same for spatial predictors
spe.rda <- dbrda(spe.hel_bray ~ netcen + updist, data = env_select)
summary(spe.rda)
##
## Call:
## dbrda(formula = spe.hel_bray ~ netcen + updist, data = env_select)
##
## Partitioning of squared Bray distance:
## Inertia Proportion
## Total 4.3288 1.0000
## Constrained 0.5154 0.1191
## Unconstrained 3.8135 0.8809
##
## Eigenvalues, and their contribution to the squared Bray distance
##
## Importance of components:
## dbRDA1 dbRDA2 MDS1 MDS2 MDS3 MDS4 MDS5 MDS6
## Eigenvalue 0.39474 0.12062 1.389 0.9772 0.7302 0.472 0.33207 0.19171
## Proportion Explained 0.09119 0.02787 0.321 0.2257 0.1687 0.109 0.07671 0.04429
## Cumulative Proportion NA NA NA NA NA NA NA NA
## MDS7 MDS8 MDS9 MDS10 MDS11 MDS12 MDS13
## Eigenvalue 0.14785 0.1195 0.09463 0.05735 0.032947 0.026192 0.015370
## Proportion Explained 0.03415 0.0276 0.02186 0.01325 0.007611 0.006051 0.003551
## Cumulative Proportion NA NA NA NA NA NA NA
## MDS14 MDS15 MDS16 iMDS1 iMDS2
## Eigenvalue 0.011064 0.0023782 0.0016478 -0.0013922 -0.0028560
## Proportion Explained 0.002556 0.0005494 0.0003807 0.0003216 0.0006598
## Cumulative Proportion NA NA NA NA NA
## iMDS3 iMDS4 iMDS5 iMDS6 iMDS7
## Eigenvalue -0.005784 -0.014384 -0.016807 -0.017753 -0.021376
## Proportion Explained 0.001336 0.003323 0.003882 0.004101 0.004938
## Cumulative Proportion NA NA NA NA NA
## iMDS8 iMDS9 iMDS10 iMDS11 iMDS12 iMDS13
## Eigenvalue -0.027152 -0.034576 -0.036679 -0.04260 -0.04684 -0.05897
## Proportion Explained 0.006272 0.007988 0.008473 0.00984 0.01082 0.01362
## Cumulative Proportion NA NA NA NA NA NA
## iMDS14 iMDS15 iMDS16 iMDS17 iMDS18
## Eigenvalue -0.06787 -0.07769 -0.08770 -0.10589 -0.12169
## Proportion Explained 0.01568 0.01795 0.02026 0.02446 0.02811
## Cumulative Proportion NA NA NA NA NA
##
## Accumulated constrained eigenvalues
## Importance of components:
## dbRDA1 dbRDA2
## Eigenvalue 0.3947 0.1206
## Proportion Explained 0.7659 0.2341
## Cumulative Proportion 0.7659 1.0000
##
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores: 3.533199
##
##
## Site scores (weighted sums of species scores)
##
## dbRDA1 dbRDA2 MDS1 MDS2 MDS3 MDS4
## 1 0.03079 0.02086 0.443572 -0.509073 -0.07769 -1.05770
## 2 0.71832 0.43116 0.474971 0.635842 -0.29326 0.50678
## 3 -0.04721 -2.74094 -0.141670 -0.088012 1.14135 -0.12995
## 4 0.01772 -3.50459 0.115171 -1.313350 0.20198 0.84532
## 5 -0.32566 -1.42176 0.654893 0.016069 0.55155 -0.06800
## 6 0.72210 -1.64050 0.570204 0.088382 0.10581 -0.10658
## 7 -1.10558 0.09717 -0.369844 -0.272616 -0.16810 0.34656
## 8 -0.18673 0.52663 0.554834 1.248341 0.06126 -0.62322
## 9 -1.87557 1.15858 -0.639086 0.781404 -0.24620 -0.14224
## 10 -0.96878 -0.76984 -0.719617 -0.290543 0.49107 0.20709
## 11 0.11261 2.48593 0.221804 0.402044 -0.80474 -0.37386
## 12 1.04844 -0.21554 0.193892 0.092443 0.25970 0.86003
## 13 -1.24808 0.44331 -0.599846 0.028805 -0.19454 0.28743
## 14 1.25611 -0.45196 1.551121 0.523263 0.10082 0.21717
## 15 -1.01413 -0.32955 -0.980964 0.093579 0.09612 -0.04228
## 16 1.90373 0.67458 0.097707 -0.008746 -0.92598 1.56757
## 17 -1.37395 0.68301 -0.612465 0.141270 -0.23712 0.34682
## 18 -1.95385 1.26346 -0.581503 0.392340 -0.24748 0.26804
## 19 0.23442 3.46276 -0.164253 0.251482 -0.83790 -0.91794
## 20 -0.12218 0.38358 -0.532349 0.230205 -0.42232 0.20971
## 21 -2.52028 1.68572 -1.325649 1.074429 0.00633 -0.22347
## 22 1.09585 -2.60626 0.251395 -0.325519 0.60016 0.15487
## 23 0.13767 -0.04685 -0.057907 -0.573964 -0.42608 -0.44829
## 24 1.81459 -0.44404 0.571144 -1.006619 -0.20588 -0.86944
## 25 -1.70361 2.75887 -0.500433 0.366802 0.32907 0.38954
## 26 1.21825 -0.62815 0.148124 0.093074 0.05832 0.64989
## 27 2.04310 1.71402 0.001172 -0.871680 -0.94463 -0.70365
## 28 -0.62007 -1.14215 -0.376631 0.105274 1.18780 -0.07127
## 29 1.34657 0.13746 1.111295 0.401562 -0.91137 0.97202
## 30 0.26768 1.28419 0.333591 0.181043 -0.17537 -0.63533
## 31 -0.45809 -2.41232 -0.198776 -0.736562 0.25061 0.32623
## 32 0.11686 -0.75156 0.344145 -0.520102 -0.03041 -0.61279
## 33 1.57525 1.84844 0.082089 -0.494682 -0.76938 0.15659
## 34 -0.47482 1.79597 -0.639963 -0.480561 -0.17134 -0.79326
## 35 0.14388 -3.64825 0.208887 -0.768867 1.43955 0.43241
## 36 0.21397 -0.19845 0.700349 1.187955 1.06427 -0.60028
## 37 -0.01932 0.09703 -0.189405 -0.074710 0.14400 -0.32451
##
##
## Site constraints (linear combinations of constraining variables)
##
## dbRDA1 dbRDA2 MDS1 MDS2 MDS3 MDS4
## 1 -0.73416 -0.17667 0.443572 -0.509073 -0.07769 -1.05770
## 2 0.02228 -0.61344 0.474971 0.635842 -0.29326 0.50678
## 3 0.79019 -0.81424 -0.141670 -0.088012 1.14135 -0.12995
## 4 -0.64513 -0.18785 0.115171 -1.313350 0.20198 0.84532
## 5 -0.70290 -0.24187 0.654893 0.016069 0.55155 -0.06800
## 6 0.33049 -0.76784 0.570204 0.088382 0.10581 -0.10658
## 7 -1.13258 0.01401 -0.369844 -0.272616 -0.16810 0.34656
## 8 0.17308 -0.73123 0.554834 1.248341 0.06126 -0.62322
## 9 0.05980 -0.66540 -0.639086 0.781404 -0.24620 -0.14224
## 10 0.11807 0.66719 -0.719617 -0.290543 0.49107 0.20709
## 11 0.18170 0.62463 0.221804 0.402044 -0.80474 -0.37386
## 12 0.52300 0.62623 0.193892 0.092443 0.25970 0.86003
## 13 -0.67811 -0.25729 -0.599846 0.028805 -0.19454 0.28743
## 14 -0.93278 -0.10672 1.551121 0.523263 0.10082 0.21717
## 15 0.44549 -0.60608 -0.980964 0.093579 0.09612 -0.04228
## 16 0.90818 -0.07702 0.097707 -0.008746 -0.92598 1.56757
## 17 -0.70842 -0.20120 -0.612465 0.141270 -0.23712 0.34682
## 18 -1.15566 0.03171 -0.581503 0.392340 -0.24748 0.26804
## 19 0.32250 0.59444 -0.164253 0.251482 -0.83790 -0.91794
## 20 0.27139 -0.57771 -0.532349 0.230205 -0.42232 0.20971
## 21 -0.02181 -0.40054 -1.325649 1.074429 0.00633 -0.22347
## 22 0.78741 -0.94715 0.251395 -0.325519 0.60016 0.15487
## 23 -0.44652 -0.38594 -0.057907 -0.573964 -0.42608 -0.44829
## 24 0.29443 -0.76841 0.571144 -1.006619 -0.20588 -0.86944
## 25 -0.49871 1.02553 -0.500433 0.366802 0.32907 0.38954
## 26 0.62527 -0.36947 0.148124 0.093074 0.05832 0.64989
## 27 1.16557 0.23795 0.001172 -0.871680 -0.94463 -0.70365
## 28 0.46140 0.52227 -0.376631 0.105274 1.18780 -0.07127
## 29 0.06443 0.75038 1.111295 0.401562 -0.91137 0.97202
## 30 -0.05978 0.76030 0.333591 0.181043 -0.17537 -0.63533
## 31 -0.45764 -0.34376 -0.198776 -0.736562 0.25061 0.32623
## 32 -0.66474 -0.22103 0.344145 -0.520102 -0.03041 -0.61279
## 33 0.59703 0.43085 0.082089 -0.494682 -0.76938 0.15659
## 34 0.11174 0.84412 -0.639963 -0.480561 -0.17134 -0.79326
## 35 0.27483 0.73037 0.208887 -0.768867 1.43955 0.43241
## 36 0.26872 0.73193 0.700349 1.187955 1.06427 -0.60028
## 37 0.04198 0.86894 -0.189405 -0.074710 0.14400 -0.32451
##
##
## Biplot scores for constraining variables
##
## dbRDA1 dbRDA2 MDS1 MDS2 MDS3 MDS4
## netcen -0.8133 0.5819 0 0 0 0
## updist -0.8636 -0.5041 0 0 0 0
anova(spe.rda)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: dbrda(formula = spe.hel_bray ~ netcen + updist, data = env_select)
## Df SumOfSqs F Pr(>F)
## Model 2 0.5154 2.2975 0.018 *
## Residual 34 3.8135
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RsquareAdj(spe.rda)$adj.r.squared
## [1] 0.06723433
mod0 <- dbrda(spe.hel_bray ~ 1, env_select[,c(9:10)]) # Model with intercept only #edit_PH
mod1 <- dbrda(spe.hel_bray ~ ., env_select[,c(9:10)]) # Model with all explanatory variables #edit_PH
step.res <- ordiR2step(mod0, mod1, direction = "both",perm.max = 200)
## Step: R2.adj= 0
## Call: spe.hel_bray ~ 1
##
## R2.adjusted
## <All variables> 0.06723433
## + updist 0.04867127
## + netcen 0.04317133
## <none> 0.00000000
##
## Df AIC F Pr(>F)
## + updist 1 54.314 2.8418 0.022 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.04867127
## Call: spe.hel_bray ~ updist
##
## R2.adjusted
## <All variables> 0.06723433
## + netcen 0.06723433
## <none> 0.04867127
##
## Df AIC F Pr(>F)
## + netcen 1 54.512 1.6965 0.172
step.res$anova # Summary table
## R2.adj Df AIC F Pr(>F)
## + updist 0.048671 1 54.314 2.8418 0.022 *
## <All variables> 0.067234
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RsquareAdj(spe.rda)$adj.r.squared
## [1] 0.06723433
anova.cca(spe.rda, step=1000);
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: dbrda(formula = spe.hel_bray ~ netcen + updist, data = env_select)
## Df SumOfSqs F Pr(>F)
## Model 2 0.5154 2.2975 0.028 *
## Residual 34 3.8135
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova.cca(spe.rda, step=1000, by="term");
## Permutation test for dbrda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## Model: dbrda(formula = spe.hel_bray ~ netcen + updist, data = env_select)
## Df SumOfSqs F Pr(>F)
## netcen 1 0.3019 2.6920 0.033 *
## updist 1 0.2134 1.9029 0.114
## Residual 34 3.8135
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RsquareAdj(spe.rda)$adj.r.squared;
## [1] 0.06723433
RsquareAdj(spe.rda)$r.squared
## [1] 0.1190546
# Generate a plot of the RDA
rownames(spe.rda$CCA$biplot) = c("network centrality","upstream distance")
ordiplot(spe.rda, choices=c(1,2), type="n", xlab="dbRDA 1 (76.59 %)", ylab="dbRDA 2 (23.41 %)")
points(spe.rda, "sites", pch=21, col="black", bg="white")
text(spe.rda, "cn", col="grey", cex=0.9, scaling=1, pos=4)
## Warning in arrows(0, 0, pts[, 1], pts[, 2], length = head.arrow, ...): "pos" is
## not a graphical parameter
## Warning in strwidth(labels, ...): "pos" is not a graphical parameter
## Warning in strheight(labels, ...): "pos" is not a graphical parameter
# Save the plot to an object
figure4d = recordPlot()
#Variation partitioning
spe.varpart1 <- varpart(spe.hel_bray, env_select[,1:8], env_select[,9:10])
plot(spe.varpart1,digits=2)
spe.varpart1
##
## Partition of squared Bray distance in dbRDA
##
## Call: varpart(Y = spe.hel_bray, X = env_select[, 1:8], env_select[,
## 9:10])
##
## Explanatory tables:
## X1: env_select[, 1:8]
## X2: env_select[, 9:10]
##
## No. of explanatory tables: 2
## Total variation (SS): 4.3288
## No. of observations: 37
##
## Partition table:
## Df R.squared Adj.R.squared Testable
## [a+c] = X1 8 0.26945 0.06073 TRUE
## [b+c] = X2 2 0.11905 0.06723 TRUE
## [a+b+c] = X1+X2 10 0.36716 0.12376 TRUE
## Individual fractions
## [a] = X1|X2 8 0.05653 TRUE
## [b] = X2|X1 2 0.06303 TRUE
## [c] 0 0.00420 FALSE
## [d] = Residuals 0.87624 FALSE
## ---
## Use function 'dbrda' to test significance of fractions of interest
anova.cca(dbrda(spe.hel_bray ~ T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + Condition(netcen + updist),
data=env_select), step=1000)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: dbrda(formula = spe.hel_bray ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander + Condition(netcen + updist), data = env_select)
## Df SumOfSqs F Pr(>F)
## Model 8 1.0740 1.2742 0.181
## Residual 26 2.7395
anova.cca(dbrda(spe.hel_bray ~ netcen + updist+
Condition(T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander), data=env_select), step=1000)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: dbrda(formula = spe.hel_bray ~ netcen + updist + Condition(T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander), data = env_select)
## Df SumOfSqs F Pr(>F)
## Model 2 0.42295 2.0071 0.069 .
## Residual 26 2.73945
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Infracommunities: Bray-Curtis dissimilarities are calculated at the individual host level Hellinger-transformed parasite data and then averaged within site
# A dummy parasite species is added to avoid problems with non-infected fishes
data_infra <- na.omit(data[,c(1,22:24,26:32)])
data_infra_disp <- dispweight(data_infra[,-1])
braycurtis <- vegdist(decostand(cbind(data_infra_disp,rep(1,nrow(data_infra))), na.rm=T, method="hellinger"), method="bray", na.rm=T)
meandist_bray <- meandist(braycurtis, data_infra[,1])
# Check whether Euclidean and Bray-Curtis distances are comparable
braycurtis <- vegdist(decostand(cbind(data_infra_disp,rep(1,nrow(data_infra))), na.rm=T, method="hellinger"), method="bray", na.rm=T)
meandist_bray <- meandist(braycurtis, data_infra[,1])
euc <- vegdist(decostand(cbind(data_infra_disp,rep(1,nrow(data_infra))), na.rm=T, method="hellinger"), method="euc", na.rm=T)
meandist_euc <- meandist(euc, data_infra[,1])
plot(meandist_bray[1:37,1:37], meandist_euc[1:37,1:37])
mantel(meandist_bray[1:37,1:37], meandist_euc[1:37,1:37])
##
## Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel(xdis = meandist_bray[1:37, 1:37], ydis = meandist_euc[1:37, 1:37])
##
## Mantel statistic r: 0.9906
## Significance: 0.001
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.177 0.234 0.297 0.346
## Permutation: free
## Number of permutations: 999
model_adonis = adonis(meandist_bray ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2)
## 'adonis' will be deprecated: use 'adonis2' instead
model_adonis$aov.tab
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## avlength 1 0.05123 0.051226 6.1160 0.11570 0.009 **
## avcondition 1 0.00393 0.003927 0.4688 0.00887 0.653
## T_av 1 0.00874 0.008741 1.0437 0.01974 0.393
## O2_sat_av 1 0.02116 0.021158 2.5262 0.04779 0.121
## Con_av 1 0.05861 0.058606 6.9971 0.13237 0.004 **
## COD_av 1 0.03543 0.035428 4.2298 0.08002 0.026 *
## NH4._av 1 0.00275 0.002754 0.3288 0.00622 0.708
## Nt_av 1 0.01049 0.010488 1.2521 0.02369 0.328
## pool_riffle 1 0.00491 0.004914 0.5867 0.01110 0.595
## meander 1 0.00959 0.009588 1.1447 0.02166 0.335
## netcen 1 0.02319 0.023187 2.7684 0.05237 0.084 .
## updist 1 0.01171 0.011707 1.3977 0.02644 0.278
## Residuals 24 0.20102 0.008376 0.45403
## Total 36 0.44274 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# environmental variables
env_select <- environment2[,c("T_av", "O2_sat_av", "Con_av", "COD_av", "NH4._av", "Nt_av", "pool_riffle", "meander", "netcen", "updist")]
env_select$pool_riffle <- as.numeric(env_select$pool_riffle)
env_select$meander <- as.numeric(env_select$meander)
pca <- prcomp(env_select, scale.=T)
summary(pca)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 1.7124 1.5545 1.1221 1.0140 0.88807 0.79463 0.56647
## Proportion of Variance 0.2933 0.2416 0.1259 0.1028 0.07887 0.06314 0.03209
## Cumulative Proportion 0.2933 0.5349 0.6608 0.7636 0.84248 0.90563 0.93771
## PC8 PC9 PC10
## Standard deviation 0.50483 0.46939 0.38429
## Proportion of Variance 0.02549 0.02203 0.01477
## Cumulative Proportion 0.96320 0.98523 1.00000
plot(pca)
biplot(pca)
# Assess the effect of environmental variables on parasite infracommunity dissimilarities using distance based RDA
spe.rda <- dbrda(meandist_bray ~ T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander, data = env_select, tidy=TRUE)
summary(spe.rda)
##
## Call:
## dbrda(formula = meandist_bray ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander, data = env_select, tidy = TRUE)
##
## Partitioning of squared Unknown distance:
## Inertia Proportion
## Total 1.5898 1.0000
## Constrained 0.4504 0.2833
## Unconstrained 1.1395 0.7167
##
## Eigenvalues, and their contribution to the squared Unknown distance
##
## Importance of components:
## dbRDA1 dbRDA2 dbRDA3 dbRDA4 dbRDA5 dbRDA6 dbRDA7
## Eigenvalue 0.2196 0.05240 0.04219 0.03708 0.02865 0.02593 0.02365
## Proportion Explained 0.1381 0.03296 0.02654 0.02333 0.01802 0.01631 0.01487
## Cumulative Proportion 0.1381 0.17107 0.19761 0.22093 0.23895 0.25527 0.27014
## dbRDA8 MDS1 MDS2 MDS3 MDS4 MDS5 MDS6
## Eigenvalue 0.02089 0.1762 0.13809 0.10605 0.06557 0.06183 0.05722
## Proportion Explained 0.01314 0.1108 0.08686 0.06671 0.04125 0.03889 0.03599
## Cumulative Proportion 0.28328 0.3941 0.48098 0.54769 0.58894 0.62783 0.66382
## MDS7 MDS8 MDS9 MDS10 MDS11 MDS12 MDS13
## Eigenvalue 0.04417 0.04047 0.03937 0.03784 0.03506 0.03412 0.03273
## Proportion Explained 0.02778 0.02545 0.02476 0.02380 0.02205 0.02146 0.02058
## Cumulative Proportion 0.69160 0.71706 0.74182 0.76562 0.78767 0.80913 0.82972
## MDS14 MDS15 MDS16 MDS17 MDS18 MDS19 MDS20
## Eigenvalue 0.02795 0.02650 0.02630 0.02513 0.02348 0.02044 0.01991
## Proportion Explained 0.01758 0.01667 0.01654 0.01580 0.01477 0.01286 0.01252
## Cumulative Proportion 0.84730 0.86397 0.88051 0.89632 0.91109 0.92394 0.93647
## MDS21 MDS22 MDS23 MDS24 MDS25 MDS26
## Eigenvalue 0.01975 0.01836 0.01714 0.015343 0.013426 0.008684
## Proportion Explained 0.01242 0.01155 0.01078 0.009651 0.008445 0.005462
## Cumulative Proportion 0.94889 0.96044 0.97122 0.980869 0.989313 0.994775
## MDS27 MDS28
## Eigenvalue 0.007381 0.0009256
## Proportion Explained 0.004642 0.0005822
## Cumulative Proportion 0.999418 1.0000000
##
## Accumulated constrained eigenvalues
## Importance of components:
## dbRDA1 dbRDA2 dbRDA3 dbRDA4 dbRDA5 dbRDA6 dbRDA7
## Eigenvalue 0.2196 0.0524 0.04219 0.03708 0.02865 0.02593 0.02365
## Proportion Explained 0.4875 0.1164 0.09368 0.08234 0.06361 0.05759 0.05251
## Cumulative Proportion 0.4875 0.6039 0.69757 0.77991 0.84353 0.90111 0.95362
## dbRDA8
## Eigenvalue 0.02089
## Proportion Explained 0.04638
## Cumulative Proportion 1.00000
##
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores: 2.750505
##
##
## Site scores (weighted sums of species scores)
##
## dbRDA1 dbRDA2 dbRDA3 dbRDA4 dbRDA5 dbRDA6
## SITE 1 -0.193401 0.063457 -0.61764 -0.312982 0.78782 -0.875368
## SITE 11 -0.424501 -0.396544 -0.49524 -0.201551 1.04886 0.509773
## SITE 12 -0.702700 0.085074 -0.27366 0.063060 0.10379 -0.004917
## SITE 13 -0.007236 -0.112556 -0.10006 0.582155 0.14248 0.722762
## SITE 14 -0.125146 0.085307 -0.29498 0.469582 -0.32057 0.196422
## SITE 15 -0.265297 -0.146273 0.10540 0.221087 -0.44493 0.501526
## SITE 16 0.017380 0.097078 -0.34338 0.027222 -0.23381 0.253728
## SITE 17 -0.143898 0.261154 -0.74786 0.425484 -1.19344 -1.632875
## SITE 18 0.694672 0.250138 -0.19428 0.116751 -0.02255 0.197381
## SITE 19 2.013583 0.700678 -0.26310 0.241350 0.28881 -0.781458
## SITE 2 -0.269160 -0.474927 -0.05506 -0.658268 -0.75242 -0.133679
## SITE 20 0.184271 -1.120553 0.05670 0.218183 -0.85772 0.307942
## SITE 21 -0.478165 0.259320 -0.31414 -0.196556 0.20691 0.132081
## SITE 22 -0.680264 0.661630 0.17236 0.300656 0.03888 0.186979
## SITE 23 0.107511 0.117471 -0.68571 0.024261 0.82469 -0.140649
## SITE 24 -0.509800 -0.585796 1.11284 0.791260 0.27686 -0.880620
## SITE 26 0.023452 0.520302 -0.75929 -0.348237 0.26634 -0.019943
## SITE 28 0.081124 0.364618 -0.32704 -0.131404 -0.21902 -0.874784
## SITE 29 0.069397 0.046767 0.74532 -1.644458 0.22297 0.183059
## SITE 3 -0.655376 -0.078746 -0.30374 0.077802 0.17383 0.171647
## SITE 30 -0.750306 0.186160 -0.53020 -0.087275 0.13772 0.020057
## SITE 31 -0.003379 -0.217568 0.16729 0.877596 0.12615 -0.182314
## SITE 32 0.525097 -0.234069 -0.10169 0.176823 -0.13431 0.273050
## SITE 33 -0.376746 -0.057152 -0.09869 0.386971 -0.79147 -0.032496
## SITE 34 -0.427146 0.315945 -0.34228 -0.588588 0.30660 -0.646267
## SITE 35 0.144043 -0.736009 1.21209 0.693519 0.92859 0.065714
## SITE 36 -0.221245 -1.014780 0.73058 -0.257720 0.04394 0.325929
## SITE 38 0.260121 0.548181 -0.12839 0.009881 0.12554 0.030680
## SITE 39 -0.912898 0.029942 -0.24674 -0.094014 0.07078 0.062174
## SITE 4 -0.039914 -0.160458 -0.52023 -1.243446 0.02703 0.155617
## SITE 40 0.437382 0.347286 -0.44241 0.551293 0.15660 0.371577
## SITE 41 1.061747 0.460820 0.41788 -0.329177 -0.20931 1.117435
## SITE 42 0.492754 -2.142555 1.38058 -0.595630 -0.70424 -0.924924
## SITE 5 0.486220 -0.097372 -0.04666 0.145692 0.78718 0.142772
## SITE 6 0.366097 0.007433 0.29450 1.020707 0.06306 1.202313
## SITE 7 -0.109776 2.157876 1.76424 -0.745567 -0.40355 -0.618902
## SITE 9 0.331503 0.008719 0.07268 0.013537 -0.86805 0.618579
##
##
## Site constraints (linear combinations of constraining variables)
##
## dbRDA1 dbRDA2 dbRDA3 dbRDA4 dbRDA5 dbRDA6
## SITE 1 -0.013706 0.06161 -0.566299 -0.56263 0.68337 -0.799399
## SITE 11 -0.498049 -0.62259 -0.387895 -0.14443 0.86286 0.337683
## SITE 12 -0.449977 -0.41292 -0.381703 0.39230 -0.03446 -0.494715
## SITE 13 0.646634 0.78752 -0.170195 0.73514 0.33780 0.894640
## SITE 14 -0.022889 0.21915 -0.559986 0.43068 -0.75852 -0.249161
## SITE 15 -0.477115 -0.18072 0.381861 -0.03891 -0.43609 0.261592
## SITE 16 -0.003779 -0.10419 -0.224779 0.04360 -0.47874 0.521790
## SITE 17 0.014080 -0.01398 -0.393448 0.37555 -0.45946 -0.612432
## SITE 18 0.433470 -0.38610 0.311804 0.49069 -0.10431 0.385204
## SITE 19 1.797552 -0.01477 -0.068343 0.14213 -0.02578 -0.648483
## SITE 2 0.097800 -0.47743 0.114446 -0.54828 -0.66048 -0.040147
## SITE 20 0.233513 -0.56819 -0.452230 0.18420 -0.43596 0.342081
## SITE 21 0.115261 0.39505 0.012004 -0.45991 0.25232 0.602002
## SITE 22 -0.271628 0.94479 0.262170 0.24479 0.06794 -0.038339
## SITE 23 -0.268975 -0.01406 -0.687107 -0.14113 0.84944 -0.150595
## SITE 24 -0.470197 -0.06626 1.074097 0.59004 0.26654 -0.627883
## SITE 26 -0.127347 0.66408 -0.992149 -0.04145 0.27143 0.182401
## SITE 28 -0.129789 0.14568 0.024292 0.34127 -0.32756 -0.565381
## SITE 29 0.193897 0.19946 0.677019 -0.93440 0.11201 0.170499
## SITE 3 -0.678081 -0.39312 -0.314876 0.16522 0.26479 0.467313
## SITE 30 -0.606829 -0.31255 -0.451988 0.29629 -0.13401 0.461075
## SITE 31 -0.139331 -0.12436 0.009754 0.63248 0.30366 -0.804481
## SITE 32 0.410116 -0.40478 -0.132639 -0.21685 -0.23921 -0.040507
## SITE 33 -0.074343 0.29536 -0.134441 0.39097 -1.04116 -0.336434
## SITE 34 -0.226204 0.29969 0.198989 -0.65741 0.28977 -0.726017
## SITE 35 0.238182 -0.36205 0.851383 0.44887 0.64367 0.065366
## SITE 36 -0.212694 -0.42956 0.325410 -0.02991 0.12247 -0.002134
## SITE 38 -0.290761 0.68825 0.099402 -0.20120 0.08529 -0.165182
## SITE 39 -0.424749 -0.34091 0.653590 -0.45291 -0.31773 0.026280
## SITE 4 0.249290 -0.18203 -0.478520 -1.12533 0.02918 0.073437
## SITE 40 0.011975 0.77178 -0.349485 0.10183 0.04876 -0.016774
## SITE 41 0.562714 0.30100 0.218979 -0.62325 -0.17919 0.463349
## SITE 42 0.137879 -0.72177 0.095094 -0.49959 -0.18024 -0.377517
## SITE 5 0.741324 -0.45934 0.254706 0.55877 0.92073 0.108541
## SITE 6 -0.173401 -0.11045 0.214506 0.34726 0.12331 0.892139
## SITE 7 -0.341411 1.02435 0.770083 -0.17660 -0.04056 -0.023937
## SITE 9 0.017570 -0.09563 0.196492 -0.05789 -0.68189 0.464125
##
##
## Biplot scores for constraining variables
##
## dbRDA1 dbRDA2 dbRDA3 dbRDA4 dbRDA5 dbRDA6
## T_av 0.19654 0.04389 0.79605 -0.13924 -0.28696 -0.21269
## O2_sat_av -0.24302 0.51209 -0.29524 -0.44557 -0.62386 -0.02169
## Con_av 0.69249 0.28619 -0.13247 -0.04714 -0.14852 -0.19236
## COD_av 0.39564 -0.11996 0.29825 0.49441 0.08302 -0.03404
## NH4._av 0.62169 -0.12508 0.08817 0.43757 0.26414 0.12573
## Nt_av 0.61125 -0.07446 -0.30978 -0.09649 0.17024 -0.41824
## pool_riffle 0.25108 0.73042 -0.13229 -0.45593 0.34283 -0.02248
## meander -0.04274 0.63037 -0.22360 0.31145 0.06195 -0.62153
rownames(spe.rda$CCA$biplot) = c("temperature","O2","conductivity","COD","NH4","Nt","pool riffle pattern","meander")
# Generate a plot of the RDA
ordiplot(spe.rda, choices=c(1,2), type="n", xlab="dbRDA 1 (48.75 %)", ylab="dbRDA 2 (11.64 %)")
points(spe.rda, "sites", pch=21, col="black", bg="white")
text(spe.rda, "cn", col="grey", cex=0.9, scaling=1)
# Save the plot to an object
figure4a = recordPlot()
anova(spe.rda)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: dbrda(formula = meandist_bray ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander, data = env_select, tidy = TRUE)
## Df SumOfSqs F Pr(>F)
## Model 8 0.45036 1.3833 0.007 **
## Residual 28 1.13946
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(spe.rda, by="term")
## Permutation test for dbrda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## Model: dbrda(formula = meandist_bray ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander, data = env_select, tidy = TRUE)
## Df SumOfSqs F Pr(>F)
## T_av 1 0.04338 1.0661 0.308
## O2_sat_av 1 0.04859 1.1940 0.237
## Con_av 1 0.12820 3.1503 0.001 ***
## COD_av 1 0.06629 1.6290 0.056 .
## NH4._av 1 0.02849 0.7002 0.835
## Nt_av 1 0.03286 0.8076 0.688
## pool_riffle 1 0.04004 0.9838 0.397
## meander 1 0.06250 1.5358 0.057 .
## Residual 28 1.13946
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RsquareAdj(spe.rda)$adj.r.squared
## [1] 0.07850103
mod0 <- dbrda(meandist_bray ~ 1, env_select) # Model with intercept only #edit_PH
mod1 <- dbrda(meandist_bray ~ ., env_select) # Model with all explanatory variables #edit_PH
step.res <- ordiR2step(mod0, mod1, direction = "both",perm.max = 200)
## Step: R2.adj= 0
## Call: meandist_bray ~ 1
##
## R2.adjusted
## <All variables> 0.0898822750
## + Con_av 0.0492715387
## + NH4._av 0.0376899022
## + Nt_av 0.0353578268
## + COD_av 0.0097867139
## + updist 0.0092890919
## + pool_riffle 0.0070398050
## + netcen 0.0034499960
## + O2_sat_av 0.0031240320
## <none> 0.0000000000
## + T_av -0.0005031246
## + meander -0.0037253892
##
## Df AIC F Pr(>F)
## + Con_av 1 17.228 2.8657 0.004 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.04927154
## Call: meandist_bray ~ Con_av
##
## R2.adjusted
## <All variables> 0.08988227
## + COD_av 0.08426096
## + NH4._av 0.07058808
## + updist 0.06800179
## + O2_sat_av 0.06266812
## + netcen 0.05455352
## + meander 0.05428744
## + Nt_av 0.05273614
## <none> 0.04927154
## + pool_riffle 0.04923534
## + T_av 0.04694754
##
## Df AIC F Pr(>F)
## + COD_av 1 16.768 2.3373 0.004 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.08426096
## Call: meandist_bray ~ Con_av + COD_av
##
## R2.adjusted
## + updist 0.09961037
## + meander 0.09248673
## <All variables> 0.08988227
## + netcen 0.08626581
## + pool_riffle 0.08559576
## <none> 0.08426096
## + T_av 0.08115910
## + O2_sat_av 0.08092067
## + Nt_av 0.07979657
## + NH4._av 0.07591794
step.res$anova # Summary table
## R2.adj Df AIC F Pr(>F)
## + Con_av 0.049272 1 17.228 2.8657 0.004 **
## + COD_av 0.084261 1 16.768 2.3373 0.004 **
## <All variables> 0.089882
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
spe.rda <- dbrda(meandist_bray ~ Con_av + COD_av, env_select)
plot(spe.rda, scaling = 1) # it is for technical reasons not possible to plot both site and species scores
summary(spe.rda)
##
## Call:
## dbrda(formula = meandist_bray ~ Con_av + COD_av, data = env_select)
##
## Partitioning of squared Unknown distance:
## Inertia Proportion
## Total 1.5898 1.0000
## Constrained 0.2148 0.1351
## Unconstrained 1.3750 0.8649
##
## Eigenvalues, and their contribution to the squared Unknown distance
##
## Importance of components:
## dbRDA1 dbRDA2 MDS1 MDS2 MDS3 MDS4 MDS5
## Eigenvalue 0.1863 0.02858 0.1972 0.15873 0.1207 0.08661 0.07076
## Proportion Explained 0.1172 0.01798 0.1240 0.09984 0.0759 0.05447 0.04451
## Cumulative Proportion NA NA NA NA NA NA NA
## MDS6 MDS7 MDS8 MDS9 MDS10 MDS11 MDS12
## Eigenvalue 0.06231 0.04802 0.04613 0.04387 0.04231 0.03826 0.03703
## Proportion Explained 0.03919 0.03021 0.02901 0.02760 0.02661 0.02406 0.02329
## Cumulative Proportion NA NA NA NA NA NA NA
## MDS13 MDS14 MDS15 MDS16 MDS17 MDS18 MDS19
## Eigenvalue 0.03557 0.03355 0.03200 0.02808 0.02797 0.02575 0.02487
## Proportion Explained 0.02237 0.02110 0.02013 0.01766 0.01759 0.01620 0.01564
## Cumulative Proportion NA NA NA NA NA NA NA
## MDS20 MDS21 MDS22 MDS23 MDS24 MDS25 MDS26
## Eigenvalue 0.02377 0.02238 0.02182 0.02104 0.01990 0.01880 0.01782
## Proportion Explained 0.01495 0.01408 0.01372 0.01323 0.01252 0.01183 0.01121
## Cumulative Proportion NA NA NA NA NA NA NA
## MDS27 MDS28 MDS29 MDS30 MDS31 MDS32
## Eigenvalue 0.01719 0.01622 0.011778 0.010499 0.007903 0.006958
## Proportion Explained 0.01081 0.01020 0.007408 0.006604 0.004971 0.004376
## Cumulative Proportion NA NA NA NA NA NA
## MDS33 iMDS1
## Eigenvalue 0.005670 -0.006426
## Proportion Explained 0.003566 0.004042
## Cumulative Proportion NA NA
##
## Accumulated constrained eigenvalues
## Importance of components:
## dbRDA1 dbRDA2
## Eigenvalue 0.1863 0.02858
## Proportion Explained 0.8670 0.13304
## Cumulative Proportion 0.8670 1.00000
##
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores: 2.750505
##
##
## Site scores (weighted sums of species scores)
##
## dbRDA1 dbRDA2 MDS1 MDS2 MDS3 MDS4
## SITE 1 -0.215744 0.44381 0.01302 0.31819 -0.28116 -0.178171
## SITE 11 -0.507248 0.82144 -0.09719 -0.28198 -0.04679 0.342258
## SITE 12 -0.746029 0.15775 0.31672 0.25334 -0.08234 -0.151368
## SITE 13 -0.005847 -0.21986 -0.02027 0.91847 -0.12233 -0.332031
## SITE 14 -0.101987 0.28502 0.01314 0.55873 -0.09066 -0.223950
## SITE 15 -0.290401 -0.51714 -0.23614 -0.10408 -0.36670 -0.263022
## SITE 16 0.016467 0.22318 0.14693 -0.09879 -0.21010 0.156923
## SITE 17 -0.053594 0.78979 0.68552 0.22623 0.55763 0.068555
## SITE 18 0.741297 -0.34140 0.34156 -0.78781 -0.33039 0.372025
## SITE 19 2.204147 0.31903 0.07576 -0.77176 -0.32619 0.237586
## SITE 2 -0.334822 0.53901 0.20181 0.23760 0.25255 0.460075
## SITE 20 0.132501 -1.43824 -0.36307 -0.15683 0.91526 0.652586
## SITE 21 -0.518164 0.15765 0.51672 0.22690 -0.05180 0.071838
## SITE 22 -0.681351 -0.71637 0.44536 0.61813 0.24285 -0.440713
## SITE 23 0.085803 -0.07384 -0.18742 -0.58764 -0.50280 0.146049
## SITE 24 -0.506368 -1.41199 -0.02580 0.53865 0.25884 0.040947
## SITE 26 0.017028 0.61108 0.37526 -0.36715 -0.10492 0.425754
## SITE 28 0.121498 0.47609 0.59944 -0.22080 -0.07107 0.465237
## SITE 29 -0.014961 0.35365 0.04914 -0.10076 -0.08378 0.641604
## SITE 3 -0.711132 -0.09645 0.08469 -0.08718 -0.21517 0.038217
## SITE 30 -0.812113 -0.12175 0.55582 -0.16918 -0.25866 0.087617
## SITE 31 0.043616 -0.61571 -0.31459 0.47663 -0.16031 -0.651535
## SITE 32 0.562390 0.35603 -0.38737 0.06837 -0.20621 -0.072343
## SITE 33 -0.372173 0.03501 0.10319 0.85337 -0.01445 -0.328696
## SITE 34 -0.466651 0.48378 0.50738 0.02349 -0.08407 0.306511
## SITE 35 0.163917 -1.66724 -0.04520 0.22652 0.57157 0.198776
## SITE 36 -0.296043 -0.76361 -0.59812 0.40106 -0.04152 -0.249681
## SITE 38 0.302702 0.34509 -0.16274 -0.49724 -0.25040 -0.005416
## SITE 39 -0.978171 0.31142 0.36167 0.40699 0.02355 0.020399
## SITE 4 -0.141634 0.44343 -0.02774 -0.11080 -0.07834 0.180628
## SITE 40 0.496377 -0.01681 -0.38350 -0.04028 -0.50240 -0.167352
## SITE 41 1.117637 0.21808 -0.82799 -0.60061 -0.49878 -0.637014
## SITE 42 0.437517 -0.97421 -1.46550 0.08641 1.33505 0.880493
## SITE 5 0.522987 0.17736 0.27513 0.18757 -0.04206 0.152978
## SITE 6 0.423826 -0.07352 -0.91744 -0.10155 -0.42928 -0.789691
## SITE 7 0.013170 1.17661 0.57227 -1.15454 1.52741 -1.616432
## SITE 9 0.351552 0.32386 -0.18044 -0.38769 -0.23204 0.160360
##
##
## Site constraints (linear combinations of constraining variables)
##
## dbRDA1 dbRDA2 MDS1 MDS2 MDS3 MDS4
## SITE 1 -0.079483 0.44275 0.01302 0.31819 -0.28116 -0.178171
## SITE 11 -0.692558 0.69501 -0.09719 -0.28198 -0.04679 0.342258
## SITE 12 -0.387516 -0.29270 0.31672 0.25334 -0.08234 -0.151368
## SITE 13 0.692033 0.28074 -0.02027 0.91847 -0.12233 -0.332031
## SITE 14 0.273072 0.58517 0.01314 0.55873 -0.09066 -0.223950
## SITE 15 -0.503973 -0.49304 -0.23614 -0.10408 -0.36670 -0.263022
## SITE 16 -0.041630 0.31054 0.14693 -0.09879 -0.21010 0.156923
## SITE 17 0.286975 0.22462 0.68552 0.22623 0.55763 0.068555
## SITE 18 0.355899 -0.83587 0.34156 -0.78781 -0.33039 0.372025
## SITE 19 1.702755 -0.40556 0.07576 -0.77176 -0.32619 0.237586
## SITE 2 -0.068404 0.54101 0.20181 0.23760 0.25255 0.460075
## SITE 20 0.005618 -0.62143 -0.36307 -0.15683 0.91526 0.652586
## SITE 21 -0.085484 -0.08616 0.51672 0.22690 -0.05180 0.071838
## SITE 22 -0.015359 -0.63175 0.44536 0.61813 0.24285 -0.440713
## SITE 23 -0.475351 -0.06131 -0.18742 -0.58764 -0.50280 0.146049
## SITE 24 -0.138656 -0.63048 -0.02580 0.53865 0.25884 0.040947
## SITE 26 -0.086478 0.70999 0.37526 -0.36715 -0.10492 0.425754
## SITE 28 0.227529 0.47771 0.59944 -0.22080 -0.07107 0.465237
## SITE 29 -0.053389 0.16940 0.04914 -0.10076 -0.08378 0.641604
## SITE 3 -0.805011 -0.44279 0.08469 -0.08718 -0.21517 0.038217
## SITE 30 -0.711357 -0.97402 0.55582 -0.16918 -0.25866 0.087617
## SITE 31 0.197349 -0.41392 -0.31459 0.47663 -0.16031 -0.651535
## SITE 32 0.342870 0.50733 -0.38737 0.06837 -0.20621 -0.072343
## SITE 33 0.260419 0.25833 0.10319 0.85337 -0.01445 -0.328696
## SITE 34 -0.213627 0.34583 0.50738 0.02349 -0.08407 0.306511
## SITE 35 0.299656 -0.71657 -0.04520 0.22652 0.57157 0.198776
## SITE 36 -0.265626 -0.39128 -0.59812 0.40106 -0.04152 -0.249681
## SITE 38 -0.177490 0.11612 -0.16274 -0.49724 -0.25040 -0.005416
## SITE 39 -0.502932 0.21901 0.36167 0.40699 0.02355 0.020399
## SITE 4 -0.213933 0.31087 -0.02774 -0.11080 -0.07834 0.180628
## SITE 40 0.157625 -0.02847 -0.38350 -0.04028 -0.50240 -0.167352
## SITE 41 0.308647 0.03805 -0.82799 -0.60061 -0.49878 -0.637014
## SITE 42 -0.043511 0.13560 -1.46550 0.08641 1.33505 0.880493
## SITE 5 0.774214 -0.17770 0.27513 0.18757 -0.04206 0.152978
## SITE 6 -0.221478 0.23059 -0.91744 -0.10155 -0.42928 -0.789691
## SITE 7 -0.080948 0.35692 0.57227 -1.15454 1.52741 -1.616432
## SITE 9 -0.020466 0.24745 -0.18044 -0.38769 -0.23204 0.160360
##
##
## Biplot scores for constraining variables
##
## dbRDA1 dbRDA2 MDS1 MDS2 MDS3 MDS4
## Con_av 0.7628 0.6467 0 0 0 0
## COD_av 0.4413 -0.8974 0 0 0 0
anova(spe.rda)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: dbrda(formula = meandist_bray ~ Con_av + COD_av, data = env_select)
## Df SumOfSqs F Pr(>F)
## Model 2 0.21484 2.6563 0.001 ***
## Residual 34 1.37498
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(spe.rda, by="term")
## Permutation test for dbrda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## Model: dbrda(formula = meandist_bray ~ Con_av + COD_av, data = env_select)
## Df SumOfSqs F Pr(>F)
## Con_av 1 0.12032 2.9752 0.001 ***
## COD_av 1 0.09452 2.3373 0.011 *
## Residual 34 1.37498
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
spe.rda <- dbrda(meandist_bray ~ netcen + updist, data = env_select)
anova(spe.rda)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: dbrda(formula = meandist_bray ~ netcen + updist, data = env_select)
## Df SumOfSqs F Pr(>F)
## Model 2 0.10742 1.2319 0.127
## Residual 34 1.48240
RsquareAdj(spe.rda)$adj.r.squared
## [1] 0.01271734
anova.cca(spe.rda, step=1000);
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: dbrda(formula = meandist_bray ~ netcen + updist, data = env_select)
## Df SumOfSqs F Pr(>F)
## Model 2 0.10742 1.2319 0.121
## Residual 34 1.48240
anova.cca(spe.rda, step=1000, by="term");
## Permutation test for dbrda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## Model: dbrda(formula = meandist_bray ~ netcen + updist, data = env_select)
## Df SumOfSqs F Pr(>F)
## netcen 1 0.04949 1.1352 0.230
## updist 1 0.05792 1.3285 0.132
## Residual 34 1.48240
RsquareAdj(spe.rda)$adj.r.squared;
## [1] 0.01271734
RsquareAdj(spe.rda)$r.squared
## [1] 0.06756637
summary(spe.rda)
##
## Call:
## dbrda(formula = meandist_bray ~ netcen + updist, data = env_select)
##
## Partitioning of squared Unknown distance:
## Inertia Proportion
## Total 1.5898 1.00000
## Constrained 0.1074 0.06757
## Unconstrained 1.4824 0.93243
##
## Eigenvalues, and their contribution to the squared Unknown distance
##
## Importance of components:
## dbRDA1 dbRDA2 MDS1 MDS2 MDS3 MDS4 MDS5
## Eigenvalue 0.05871 0.04871 0.3207 0.1666 0.10533 0.08746 0.06478
## Proportion Explained 0.03693 0.03064 0.2017 0.1048 0.06625 0.05501 0.04075
## Cumulative Proportion NA NA NA NA NA NA NA
## MDS6 MDS7 MDS8 MDS9 MDS10 MDS11 MDS12
## Eigenvalue 0.06249 0.04788 0.04637 0.04324 0.04235 0.03756 0.03592
## Proportion Explained 0.03931 0.03012 0.02917 0.02720 0.02664 0.02362 0.02260
## Cumulative Proportion NA NA NA NA NA NA NA
## MDS13 MDS14 MDS15 MDS16 MDS17 MDS18 MDS19
## Eigenvalue 0.03460 0.03337 0.03217 0.02941 0.02694 0.02620 0.02508
## Proportion Explained 0.02176 0.02099 0.02023 0.01850 0.01695 0.01648 0.01577
## Cumulative Proportion NA NA NA NA NA NA NA
## MDS20 MDS21 MDS22 MDS23 MDS24 MDS25 MDS26
## Eigenvalue 0.02379 0.02246 0.02179 0.02106 0.02039 0.01971 0.01801
## Proportion Explained 0.01496 0.01413 0.01371 0.01325 0.01282 0.01240 0.01133
## Cumulative Proportion NA NA NA NA NA NA NA
## MDS27 MDS28 MDS29 MDS30 MDS31 MDS32
## Eigenvalue 0.01620 0.013849 0.010275 0.009400 0.008019 0.007311
## Proportion Explained 0.01019 0.008711 0.006463 0.005913 0.005044 0.004599
## Cumulative Proportion NA NA NA NA NA NA
## MDS33 iMDS1
## Eigenvalue 0.006389 -0.004652
## Proportion Explained 0.004019 0.002926
## Cumulative Proportion NA NA
##
## Accumulated constrained eigenvalues
## Importance of components:
## dbRDA1 dbRDA2
## Eigenvalue 0.05871 0.04871
## Proportion Explained 0.54653 0.45347
## Cumulative Proportion 0.54653 1.00000
##
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores: 2.750505
##
##
## Site scores (weighted sums of species scores)
##
## dbRDA1 dbRDA2 MDS1 MDS2 MDS3 MDS4
## SITE 1 -0.59427 -0.16451 -0.14771 -0.33607 -0.0859223 -0.212994
## SITE 11 -0.45337 -0.46011 -0.23978 0.15386 0.1261509 0.237370
## SITE 12 -0.54652 -0.75642 -0.44877 0.23840 -0.2252134 -0.115066
## SITE 13 -0.12157 -0.51390 0.05860 -0.62753 -0.0424431 -0.256505
## SITE 14 -0.31417 -0.28591 -0.04140 -0.38730 0.0396086 -0.203654
## SITE 15 -0.30246 -0.79945 -0.03439 -0.11317 -0.1309347 -0.211263
## SITE 16 -0.45766 0.40537 -0.02829 -0.07163 0.0308933 0.125735
## SITE 17 -0.50007 -0.43325 -0.06301 0.73800 0.4760396 0.361411
## SITE 18 -0.23591 1.09709 0.61113 0.80322 -0.2159329 0.452248
## SITE 19 0.55697 2.59188 1.51522 0.72623 -0.4623505 0.265834
## SITE 2 0.07451 0.03855 -0.38949 0.04156 -0.0701938 0.374170
## SITE 20 1.27192 0.05240 0.05244 -0.18703 0.8794177 0.804141
## SITE 21 -0.64222 -0.20393 -0.40958 0.09451 -0.0002606 0.117485
## SITE 22 -0.69477 -0.39862 -0.55866 -0.39024 0.5266838 -0.236789
## SITE 23 -0.23341 -0.07011 0.19956 0.31332 -0.3893083 0.086348
## SITE 24 0.41401 -1.20918 -0.35479 -0.23012 0.0081164 0.193247
## SITE 26 -0.53218 0.71759 -0.02716 0.45515 -0.0183075 0.359351
## SITE 28 -0.61852 1.00385 -0.02689 0.41929 0.0029756 0.430581
## SITE 29 0.47626 0.43365 -0.12521 0.12478 -0.5950115 0.342803
## SITE 3 -0.54249 -0.73219 -0.44976 0.05668 -0.1304354 0.046421
## SITE 30 -0.88529 -0.46309 -0.61315 0.33794 -0.2093706 0.150316
## SITE 31 0.15702 -0.98675 0.23967 -0.28872 -0.1116307 -0.524180
## SITE 32 0.06463 0.23715 0.51475 -0.29819 0.0294093 -0.097103
## SITE 33 -0.24919 -0.97106 -0.16670 -0.25156 -0.1300399 -0.272378
## SITE 34 -0.53685 0.34384 -0.58079 0.14271 -0.2847025 0.213550
## SITE 35 0.57845 -0.61546 0.17558 0.01290 0.5881978 0.538055
## SITE 36 1.13941 -1.04701 -0.16168 -0.55971 -0.5621101 -0.496346
## SITE 38 0.11771 0.53367 0.16337 0.22829 -0.3411253 -0.095943
## SITE 39 -0.71027 -0.83212 -0.87000 -0.10688 -0.2047559 -0.030265
## SITE 4 0.14585 0.39575 -0.22209 -0.06629 -0.2564930 0.011573
## SITE 40 -0.18800 0.14387 0.46789 -0.34817 -0.1578324 -0.142457
## SITE 41 0.07461 0.82427 1.00639 -0.38424 0.1709964 -0.674233
## SITE 42 2.33398 -0.53599 0.44914 -1.03948 1.4092795 0.827914
## SITE 5 0.27292 0.89196 0.20209 0.11649 -0.4067229 0.125238
## SITE 6 0.45087 -0.25704 0.32728 -0.67907 -0.2681291 -0.876440
## SITE 7 0.87161 1.21367 -0.16936 1.32868 1.3928171 -1.620007
## SITE 9 0.35845 0.81153 0.14558 0.03335 -0.3813597 0.001832
##
##
## Site constraints (linear combinations of constraining variables)
##
## dbRDA1 dbRDA2 MDS1 MDS2 MDS3 MDS4
## SITE 1 -0.53436 0.244983 -0.14771 -0.33607 -0.0859223 -0.212994
## SITE 11 -0.28142 -0.386202 -0.23978 0.15386 0.1261509 0.237370
## SITE 12 0.09204 -0.878470 -0.44877 0.23840 -0.2252134 -0.115066
## SITE 13 -0.48524 0.195319 0.05860 -0.62753 -0.0424431 -0.256505
## SITE 14 -0.54659 0.190041 -0.04140 -0.38730 0.0396086 -0.203654
## SITE 15 -0.16703 -0.628956 -0.03439 -0.11317 -0.1309347 -0.211263
## SITE 16 -0.68651 0.553341 -0.02829 -0.07163 0.0308933 0.125735
## SITE 17 -0.24577 -0.530833 -0.06301 0.73800 0.4760396 0.361411
## SITE 18 -0.28345 -0.436057 0.61113 0.80322 -0.2159329 0.452248
## SITE 19 0.39319 0.351596 1.51522 0.72623 -0.4623505 0.265834
## SITE 2 0.41166 0.294943 -0.38949 0.04156 -0.0701938 0.374170
## SITE 20 0.62134 0.131756 0.05244 -0.18703 0.8794177 0.804141
## SITE 21 -0.53883 0.168679 -0.40958 0.09451 -0.0002606 0.117485
## SITE 22 -0.62229 0.383342 -0.55866 -0.39024 0.5266838 -0.236789
## SITE 23 -0.01883 -0.585259 0.19956 0.31332 -0.3893083 0.086348
## SITE 24 0.51885 -0.483973 -0.35479 -0.23012 0.0081164 0.193247
## SITE 26 -0.53040 0.217589 -0.02716 0.45515 -0.0183075 0.359351
## SITE 28 -0.69213 0.575274 -0.02689 0.41929 0.0029756 0.430581
## SITE 29 0.48333 0.208739 -0.12521 0.12478 -0.5950115 0.342803
## SITE 3 -0.11176 -0.484151 -0.44976 0.05668 -0.1304354 0.046421
## SITE 30 -0.20601 -0.234679 -0.61315 0.33794 -0.2093706 0.150316
## SITE 31 0.02640 -0.958494 0.23967 -0.28872 -0.1116307 -0.524180
## SITE 32 -0.45895 -0.021464 0.51475 -0.29819 0.0294093 -0.097103
## SITE 33 -0.18938 -0.611963 -0.16670 -0.25156 -0.1300399 -0.272378
## SITE 34 0.18801 0.867602 -0.58079 0.14271 -0.2847025 0.213550
## SITE 35 0.20501 -0.526901 0.17558 0.01290 0.5881978 0.538055
## SITE 36 0.82790 -0.414978 -0.16168 -0.55971 -0.5621101 -0.496346
## SITE 38 0.53363 0.097754 0.16337 0.22829 -0.3411253 -0.095943
## SITE 39 0.40037 0.428319 -0.87000 -0.10688 -0.2047559 -0.030265
## SITE 4 0.32911 0.494131 -0.22209 -0.06629 -0.2564930 0.011573
## SITE 40 -0.44547 0.009707 0.46789 -0.34817 -0.1578324 -0.142457
## SITE 41 -0.51320 0.184444 1.00639 -0.38424 0.1709964 -0.674233
## SITE 42 0.57268 -0.023441 0.44914 -1.03948 1.4092795 0.827914
## SITE 5 0.47441 0.462943 0.20209 0.11649 -0.4067229 0.125238
## SITE 6 0.51953 0.314869 0.32728 -0.67907 -0.2681291 -0.876440
## SITE 7 0.51653 0.318764 -0.16936 1.32868 1.3928171 -1.620007
## SITE 9 0.44365 0.511687 0.14558 0.03335 -0.3813597 0.001832
##
##
## Biplot scores for constraining variables
##
## dbRDA1 dbRDA2 MDS1 MDS2 MDS3 MDS4
## netcen -0.2800 0.9600 0 0 0 0
## updist -0.9905 0.1372 0 0 0 0
# Generate a plot of the RDA
rownames(spe.rda$CCA$biplot) = c("network centrality","upstream distance")
ordiplot(spe.rda, choices=c(1,2), type="n", xlab="dbRDA 1 (54.65 %)", ylab="dbRDA 2 (45.35 %)")
points(spe.rda, "sites", pch=21, col="black", bg="white")
text(spe.rda, "cn", col="grey", cex=0.9, scaling=1, pos=4)
## Warning in arrows(0, 0, pts[, 1], pts[, 2], length = head.arrow, ...): "pos" is
## not a graphical parameter
## Warning in strwidth(labels, ...): "pos" is not a graphical parameter
## Warning in strheight(labels, ...): "pos" is not a graphical parameter
# Save the plot to an object
figure4b = recordPlot()
#Variation partitioning
spe.varpart1 <- varpart(meandist_bray, env_select[,1:8], env_select[,9:10])
plot(spe.varpart1,digits=2)
spe.varpart1
##
## Partition of squared Unknown user-supplied distance in dbRDA
##
## Call: varpart(Y = meandist_bray, X = env_select[, 1:8], env_select[,
## 9:10])
##
## Explanatory tables:
## X1: env_select[, 1:8]
## X2: env_select[, 9:10]
##
## No. of explanatory tables: 2
## Total variation (SS): 1.5898
## No. of observations: 37
##
## Partition table:
## Df R.squared Adj.R.squared Testable
## [a+c] = X1 8 0.28328 0.07850 TRUE
## [b+c] = X2 2 0.06757 0.01272 TRUE
## [a+b+c] = X1+X2 10 0.34269 0.08988 TRUE
## Individual fractions
## [a] = X1|X2 8 0.07716 TRUE
## [b] = X2|X1 2 0.01138 TRUE
## [c] 0 0.00134 FALSE
## [d] = Residuals 0.91012 FALSE
## ---
## Use function 'dbrda' to test significance of fractions of interest
anova.cca(dbrda(meandist_bray ~ T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + Condition(netcen + updist),
data=env_select), step=1000)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: dbrda(formula = meandist_bray ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander + Condition(netcen + updist), data = env_select)
## Df SumOfSqs F Pr(>F)
## Model 8 0.4374 1.3603 0.013 *
## Residual 26 1.0450
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova.cca(dbrda(meandist_bray ~ netcen + updist+
Condition(T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander), data=env_select), step=1000)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: dbrda(formula = meandist_bray ~ netcen + updist + Condition(T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander), data = env_select)
## Df SumOfSqs F Pr(>F)
## Model 2 0.09446 1.1751 0.176
## Residual 26 1.04500
library(gridExtra)
#pdf("Figure4.pdf", width = 10, height = 10)
#grid.arrange(figure4a, figure4b, figure4c, figure4d)
#dev.off()
library(ggpubr)
library(gridGraphics)
## Loading required package: grid
pdf("Figure4.pdf", width = 16, height = 10)
ggarrange(plotlist = list(figure4a, figure4c, figure4b, figure4d), nrow = 2, ncol=2)
dev.off()
## png
## 2
# Conduct a PCA on the six physico-chemical variables from VMM
PCA = prcomp(environment2[, c(2:7)], center = TRUE, scale. = TRUE)
summary(PCA) # Summary statistics of the PCA
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6
## Standard deviation 1.5652 1.2707 0.9683 0.66323 0.62754 0.40523
## Proportion of Variance 0.4083 0.2691 0.1563 0.07331 0.06563 0.02737
## Cumulative Proportion 0.4083 0.6774 0.8337 0.90700 0.97263 1.00000
biplot(PCA) # Graphical representation of the PCA
var <- get_pca_var(PCA)
var$contrib # Show contribution of variables to PCs
## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5 Dim.6
## T_av 4.781229 3.2285483 86.3186845 1.24296167 3.7819969 0.6465796
## O2_sat_av 21.325906 12.6315198 4.5711032 37.97350069 10.9404721 12.5574980
## Con_av 1.259548 46.2188894 1.3713634 23.59374367 26.5347837 1.0216720
## COD_av 20.367872 15.5236356 1.6185615 9.07873660 46.8437587 6.5674354
## NH4._av 35.988850 0.3815385 0.8826774 0.05274591 0.3244479 62.3697404
## Nt_av 16.276595 22.0158684 5.2376100 28.05831146 11.5745407 16.8370745
pc1 = PCA$x[,1]
data_p = data
levels(data_p$site) = as.factor(environment2$Site)
environment3 = cbind(environment2, pc1)
m = merge(x = environment3, y = data_p, by.x = "Site", by.y = "site", all.y = TRUE)
m2 = m[complete.cases(m[, c("SMI", "PI", "PI_ecto", "PI_endo")]), ]
model_PI = psem(
lme(PI ~ pc1 + SMI, random = ~1|Site, data=m2),
lme(SMI ~ pc1, random = ~1|Site, data=m2)
)
plot(model_PI)
summary(model_PI)
##
## Structural Equation Model of model_PI
##
## Call:
## PI ~ pc1 + SMI
## SMI ~ pc1
##
## AIC BIC
## 18.000 57.079
##
## ---
## Tests of directed separation:
##
## No independence claims present. Tests of directed separation not possible.
##
## Global goodness-of-fit:
##
## Fisher's C = 0 with P-value = 1 and on 0 degrees of freedom
##
## ---
## Coefficients:
##
## Response Predictor Estimate Std.Error DF Crit.Value P.Value Std.Estimate
## PI pc1 0.2698 0.1154 35 2.3376 0.0253 0.1687 *
## PI SMI 0.1889 0.3794 530 0.4980 0.6187 0.0199
## SMI pc1 0.0040 0.0088 35 0.4590 0.6491 0.0239
##
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05
##
## ---
## Individual R-squared:
##
## Response method Marginal Conditional
## PI none 0.03 0.17
## SMI none 0.00 0.04
model_PI_ecto = psem(
lme(PI_ecto ~ pc1 + SMI, random = ~1|Site, data=m2),
lme(SMI ~ pc1, random = ~1|Site, data=m2)
)
plot(model_PI_ecto)
summary(model_PI_ecto)
##
## Structural Equation Model of model_PI_ecto
##
## Call:
## PI_ecto ~ pc1 + SMI
## SMI ~ pc1
##
## AIC BIC
## 18.000 57.079
##
## ---
## Tests of directed separation:
##
## No independence claims present. Tests of directed separation not possible.
##
## Global goodness-of-fit:
##
## Fisher's C = 0 with P-value = 1 and on 0 degrees of freedom
##
## ---
## Coefficients:
##
## Response Predictor Estimate Std.Error DF Crit.Value P.Value Std.Estimate
## PI_ecto pc1 0.3092 0.0851 35 3.6342 0.0009 0.2812 ***
## PI_ecto SMI 0.1335 0.2525 530 0.5286 0.5973 0.0205
## SMI pc1 0.0040 0.0088 35 0.4590 0.6491 0.0239
##
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05
##
## ---
## Individual R-squared:
##
## Response method Marginal Conditional
## PI_ecto none 0.08 0.25
## SMI none 0.00 0.04
model_PI_endo = psem(
lme(PI_endo ~ pc1 + SMI, random = ~1|Site, data=m2),
lme(SMI ~ pc1, random = ~1|Site, data=m2)
)
plot(model_PI_endo)
summary(model_PI_endo)
##
## Structural Equation Model of model_PI_endo
##
## Call:
## PI_endo ~ pc1 + SMI
## SMI ~ pc1
##
## AIC BIC
## 18.000 57.079
##
## ---
## Tests of directed separation:
##
## No independence claims present. Tests of directed separation not possible.
##
## Global goodness-of-fit:
##
## Fisher's C = 0 with P-value = 1 and on 0 degrees of freedom
##
## ---
## Coefficients:
##
## Response Predictor Estimate Std.Error DF Crit.Value P.Value Std.Estimate
## PI_endo pc1 -0.0538 0.0792 35 -0.6792 0.5015 -0.0454
## PI_endo SMI 0.1217 0.2894 530 0.4205 0.6743 0.0173
## SMI pc1 0.0040 0.0088 35 0.4590 0.6491 0.0239
##
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05
##
## ---
## Individual R-squared:
##
## Response method Marginal Conditional
## PI_endo none 0 0.11
## SMI none 0 0.04